Ecology

Moritz, M. A. et al. Learning to coexist with wildfire. Nature 515, 58–66 (2014).ADS 
CAS 
PubMed 

Google Scholar 
Bowman, D. M. et al. Vegetation fires in the Anthropocene. Nat. Rev. Earth Environ. 1, 500–515 (2020).ADS 

Google Scholar 
Fischer, A. P. et al. Wildfire risk as a socioecological pathology. Front. Ecol. Environ. 14, 276–284 (2016).
Google Scholar 
Steffensen, V. Fire Country: How Indigenous Fire Management Could Help Save Australia (Hardie Grant Books, 2020).
Google Scholar 
Australian Government, Royal Commission into National Natural Disaster Arrangements. Commonwealth Letters Patent—20 February, 2020. https://naturaldisaster.royalcommission.gov.au/publications/commonwealth-letters-patent-20-february-2020 (2020).Bowman, D. M. The impact of Aboriginal landscape burning on the Australian biota. New Phytol. 140, 385–410 (1998).CAS 
PubMed 

Google Scholar 
Roos, C. I., Williamson, G. J. & Bowman, D. M. Is anthropogenic pyrodiversity invisible in paleofire records?. Fire 2, 42 (2019).
Google Scholar 
Liebmann, M. J. et al. Native American depopulation, reforestation, and fire regimes in the Southwest United States, 1492–1900 CE. Proc. Natl. Acad. Sci. 113, E696–E704 (2016).CAS 
PubMed 
PubMed Central 

Google Scholar 
Fletcher, M.-S., Hamilton, R., Dressler, W. & Palmer, L. Indigenous knowledge and the shackles of wilderness. Proc. Natl. Acad. Sci. 118, e2022218118 (2021).CAS 
PubMed 
PubMed Central 

Google Scholar 
Thomson, D. F. Arnhem land: Explorations among an unknown people part I. The journey to Bennet Bay. Geogr. J. 112, 146–164. https://doi.org/10.2307/1789695 (1948).Article 

Google Scholar 
Yibarbuk, D. et al. Fire ecology and Aboriginal land management in central Arnhem Land, northern Australia: A tradition of ecosystem management. J. Biogeogr. 28, 325–343 (2001).
Google Scholar 
Jones, G. M. & Tingley, M. W. Pyrodiversity and biodiversity: A history, synthesis, and outlook. Divers. Distrib. 28, 386–403 (2021).
Google Scholar 
Steel, Z. L., Collins, B. M., Sapsis, D. B. & Stephens, S. L. Quantifying pyrodiversity and its drivers. Proc. R. Soc. B 288, 20203202 (2021).PubMed 
PubMed Central 

Google Scholar 
Bird, R. B., Tayor, N., Codding, B. F. & Bird, D. W. Niche construction and Dreaming logic: Aboriginal patch mosaic burning and varanid lizards (Varanus gouldii) in Australia. Proc. R. Soc. B: Biol. Sci. 280, 20132297 (2013).
Google Scholar 
Bowman, D. M., Walsh, A. & Prior, L. D. Landscape analysis of Aboriginal fire management in Central Arnhem Land, north Australia. J. Biogeogr. 31, 207–223 (2004).
Google Scholar 
Haynes, C. D. in Proceedings of the Ecological Society of Australia (Australia) (Darwin Institute of Technology).Murphy, B. P. & Bowman, D. M. The interdependence of fire, grass, kangaroos and Australian Aborigines: A case study from central Arnhem Land, northern Australia. J. Biogeogr. 34, 237–250 (2007).
Google Scholar 
Bowman, D., Garde, M. & Saulwick, A. in Histories of Old Ages: Essays in Honour of Rhys Jones (eds Anderson, A. et al.) 61–78 (Australian National University, 2001).Bowman, D. & Panton, W. Decline of Callitris intratropica RT Baker & HG Smith in the Northern Territory: Implications for pre-and post-European colonization fire regimes. J. Biogeogr. 20, 373–381 (1993).
Google Scholar 
Sharp, B. R. & Bowman, D. M. Patterns of long-term woody vegetation change in a sandstone-plateau savanna woodland, Northern Territory, Australia. J. Trop. Ecol. 20, 259–270 (2004).
Google Scholar 
Trauernicht, C., Murphy, B. P., Tangalin, N. & Bowman, D. M. Cultural legacies, fire ecology, and environmental change in the Stone Country of Arnhem Land and Kakadu National Park, Australia. Ecol. Evol. 3, 286–297 (2013).PubMed 
PubMed Central 

Google Scholar 
Edwards, A. C. & Russell-Smith, J. Ecological thresholds and the status of fire-sensitive vegetation in western Arnhem Land, northern Australia: Implications for management. Int. J. Wildland Fire 18, 127–146 (2009).
Google Scholar 
Yates, C. & Russell-Smith, J. Fire regimes and vegetation sensitivity analysis: An example from Bradshaw Station, monsoonal northern Australia. Int. J. Wildland Fire 12, 349–358 (2003).
Google Scholar 
McVicar, D. Reports Concerning Marketable Timbers and Forest Products of Several Regions of the North-West Part of the State (WA Forests Department, 1922).
Google Scholar 
Bowman, D. M., Price, O., Whitehead, P. J. & Walsh, A. The ‘wilderness effect’ and the decline of Callitris intratropica on the Arnhem Land Plateau, northern Australia. Aust. J. Bot. 49, 665–672 (2001).
Google Scholar 
Prior, L. D., McCaw, W. L., Grierson, P. F., Murphy, B. P. & Bowman, D. M. Population structures of the widespread Australian conifer Callitris columellaris are a bio-indicator of continental environmental change. For. Ecol. Manag. 262, 252–262 (2011).
Google Scholar 
Bowman, D. M., MacDermott, H. J., Nichols, S. C. & Murphy, B. P. A grass–fire cycle eliminates an obligate-seeding tree in a tropical savanna. Ecol. Evol. 4, 4185–4194 (2014).PubMed 
PubMed Central 

Google Scholar 
Lawes, M. J., Richards, A., Dathe, J. & Midgley, J. J. Bark thickness determines fire resistance of selected tree species from fire-prone tropical savanna in north Australia. Plant Ecol. 212, 2057–2069 (2011).
Google Scholar 
Bowman, D. M., Haverkamp, C., Rann, K. D. & Prior, L. D. Differential demographic filtering by surface fires: How fuel type and fuel load affect sapling mortality of an obligate seeder savanna tree. J. Ecol. 106, 1010–1022 (2018).
Google Scholar 
Trauernicht, C., Murphy, B. P., Prior, L. D., Lawes, M. J. & Bowman, D. M. Human-imposed, fine-grained patch burning explains the population stability of a fire-sensitive conifer in a frequently burnt northern Australia savanna. Ecosystems 19, 896–909 (2016).
Google Scholar 
Bininj Kunwok Regional Language Centre. https://bininjkunwok.org.au (2021).Cooke, P. M. Buffalo and tin, baki and Jesus. In Culture, Ecology and Economy of Fire Management in North Australian Savannas: Rekindling the Wurrk Tradition (eds Russell-Smith, J. et al.) 69–83 (Csiro Publishing, 2009).
Google Scholar 
Edwards, A. et al. Transforming fire management in northern Australia through successful implementation of savanna burning emissions reductions projects. J. Environ. Manag. 290, 112568 (2021).
Google Scholar 
Clarkson, C. et al. Human occupation of northern Australia by 65,000 years ago. Nature 547, 306–310 (2017).ADS 
CAS 
PubMed 

Google Scholar 
Tobler, R. et al. Aboriginal mitogenomes reveal 50,000 years of regionalism in Australia. Nature 544, 180–184 (2017).ADS 
CAS 
PubMed 

Google Scholar 
Press, T., Lea, D., Webb, A. & Alistair, G. Kakadu Natural and Cultural Heritage and Management (The Australian National University, 1995).
Google Scholar 
Murphy, B. P., Cochrane, M. A. & Russell-Smith, J. Prescribed burning protects endangered tropical heathlands of the Arnhem Plateau, northern Australia. J. Appl. Ecol. 52, 980–991 (2015).
Google Scholar 
Russell-Smith, J., Needham, S. & Brock, J. in Kakadu: Natural and Cultural Heritage and Management (eds A. J. Press et al.) 128–166 (Australian National University, 1995).Haynes, C., Ridpath, M. & Williams, M. A. Monsoonal Australia: Landscape, Ecology and Man in Northern Lowlands (CRC Press, 1991).
Google Scholar 
Bowman, D. M. et al. Biogeography of the Australian monsoon tropics. J. Biogeogr. 37, 201–216 (2010).
Google Scholar 
Williamson, G. J. et al. Measurement of inter-and intra-annual variability of landscape fire activity at a continental scale: The Australian case. Environ. Res. Lett. 11, 035003 (2016).ADS 

Google Scholar 
Russell-Smith, J. et al. Bushfires down under: Patterns and implications of contemporary Australian landscape burning. Int. J. Wildland Fire 16, 361–377. https://doi.org/10.1071/WF07018 (2007).Article 

Google Scholar 
Evans, J. & Russell-Smith, J. Delivering effective savanna fire management for defined biodiversity conservation outcomes: An Arnhem Land case study. Int. J. Wildland Fire 29, 386–400 (2019).
Google Scholar 
Corey, B. et al. Better biodiversity accounting is needed to prevent bioperversity and maximize co-benefits from savanna burning. Conserv. Lett. 13, e12685 (2020).
Google Scholar 
Crisp, M. D. et al. Turnover of southern cypresses in the post-Gondwanan world: Extinction, transoceanic dispersal, adaptation and rediversification. New Phytol. 221, 2308–2319 (2019).PubMed 

Google Scholar 
Prior, L. D. & Bowman, D. M. Classification of post-fire responses of woody plants to include pyrophobic communities. Fire 3, 15 (2020).
Google Scholar 
Brodribb, T. J. et al. Conservative water management in the widespread conifer genus Callitris. AoB Plants 5, plt052 (2013).PubMed Central 

Google Scholar 
Sakaguchi, S. et al. Climate, not Aboriginal landscape burning, controlled the historical demography and distribution of fire-sensitive conifer populations across Australia. Proc. R. Soc. B: Biol. Sci. 280, 20132182 (2013).
Google Scholar 
Allen, K. J. et al. Two climate-sensitive tree-ring chronologies from Arnhem Land, monsoonal Australia. Austral Ecol. 44, 581–596 (2019).
Google Scholar 
Baker, P. J., Palmer, J. G. & D’Arrigo, R. The dendrochronology of Callitris intratropica in northern Australia: Annual ring structure, chronology development and climate correlations. Aust. J. Bot. 56, 311–320 (2008).
Google Scholar 
Hammer, G. Site classification and tree diameter-height-age relationships for cypress pine in the Top End of the Northern Territory. Aust. For. 44, 35–41 (1981).ADS 

Google Scholar 
Prior, L., Bowman, D. & Brook, B. Growth and survival of two north Australian relictual tree species, Allosyncarpia ternata (Myrtaceae) and Callitris intratropica (Cupressaceae). Ecol. Res. 22, 228–236 (2007).
Google Scholar 
Stocker, G. Aspects of the Seeding Habits of Callitris intratropica (Forestry and Timber Bureau, 1966).
Google Scholar 
Bowman, D., Wilson, B. & Davis, G. Response of Callitris intratropica RT Baker & HG Smith to fire protection, Murgenella, northern Australia. Aust. J. Ecol. 13, 147–159 (1988).
Google Scholar 
Hawkins, P. Seed production and litter fall studies of Callitris columellaris. Aust. For. Res. 2, 3–16 (1966).
Google Scholar 
Lawes, M. J., Taplin, P., Bellairs, S. M. & Franklin, D. C. A trade-off in stand size effects in the reproductive biology of a declining tropical conifer Callitris intratropica. Plant Ecol. 214, 169–174 (2013).
Google Scholar 
Petty, A. M. & Bowman, D. M. A satellite analysis of contrasting fire patterns in aboriginal-and euro-Australian lands in tropical North Australia. Fire Ecol. 3, 32–47 (2007).
Google Scholar 
Bowman, D. & Prior, L. Impact of Aboriginal landscape burning on woody vegetation in Eucalyptus tetrodonta savanna in Arnhem Land, northern Australia. J. Biogeogr. 31, 807–817 (2004).
Google Scholar 
Trauernicht, C., Murphy, B. P., Portner, T. E. & Bowman, D. M. Tree cover–fire interactions promote the persistence of a fire-sensitive conifer in a highly flammable savanna. J. Ecol. 100, 958–968 (2012).
Google Scholar 
Russell-Smith, J. Recruitment dynamics of the long-lived obligate seeders Callitris intratropica (Cupressaceae) and Petraeomyrtus punicea (Myrtaceae). Aust. J. Bot. 54, 479–485 (2006).
Google Scholar 
Trauernicht, C., Brook, B. W., Murphy, B. P., Williamson, G. J. & Bowman, D. M. Local and global pyrogeographic evidence that indigenous fire management creates pyrodiversity. Ecol. Evol. 5, 1908–1918 (2015).PubMed 
PubMed Central 

Google Scholar 
D’Antonio, C. M. & Vitousek, P. M. Biological invasions by exotic grasses, the grass/fire cycle, and global change. Annu. Rev. Ecol. Syst. 23, 63–87 (1992).
Google Scholar 
Bowman, D. M., Franklin, D. C., Price, O. F. & Brook, B. W. Land management affects grass biomass in the Eucalyptus tetrodonta savannas of monsoonal Australia. Austral Ecol. 32, 446–452 (2007).
Google Scholar 
Cochrane, M. A. & Bowman, D. M. Manage fire regimes, not fires. Nat. Geosci. 14, 1–3 (2021).
Google Scholar 
Huffman, M. R. The many elements of traditional fire knowledge: Synthesis, classification, and aids to cross-cultural problem solving in fire-dependent systems around the world. Ecol. Soc. 18, 3 (2013).
Google Scholar 
Roos, C. I. et al. Native American fire management at an ancient wildland–urban interface in the Southwest United States. Proc. Natl. Acad. Sci. 118, e2018733118 (2021).CAS 
PubMed 
PubMed Central 

Google Scholar 
Long, J. W., Lake, F. K. & Goode, R. W. The importance of Indigenous cultural burning in forested regions of the Pacific West, USA. For. Ecol. Manag. 500, 119597 (2021).
Google Scholar 
Lake, F. K. et al. Returning fire to the land: Celebrating traditional knowledge and fire. J. For. 115, 343–353 (2017).ADS 

Google Scholar 
Petty, A. M., deKoninck, V. & Orlove, B. Cleaning, protecting, or abating? Making indigenous fire management “work” in northern Australia. J. Ethnobiol. 35, 140–162 (2015).
Google Scholar 
Bird, R. B. & Nimmo, D. Restore the lost ecological functions of people. Nat. Ecol. Evol. 2, 1050–1052 (2018).
Google Scholar 
Bowman, D. M. & Legge, S. Pyrodiversity—Why managing fire in food webs is relevant to restoration ecology. Restor. Ecol. 24, 848–853 (2016).
Google Scholar 
Trisos, C. H., Auerbach, J. & Katti, M. Decoloniality and anti-oppressive practices for a more ethical ecology. Nat. Ecol. Evol. 5, 1–8 (2021).
Google Scholar 
Department of Environment and Science, Q. G. Seasonal Surface Reflectance—Landsat, JRSRP Algorithm, Australia Coverage Dataset Version 1.0.0. https://portal.tern.org.au/seasonal-surface-reflectance-australia-coverage/22021 (2014).Key, C. & Benson, N. Landscape Assessment (LA) Sampling and Analysis Methods (USDA Forest Service, 2006).
Google Scholar 
Edwards, A. C., Russell-Smith, J. & Maier, S. W. A comparison and validation of satellite-derived fire severity mapping techniques in fire prone north Australian savannas: Extreme fires and tree stem mortality. Remote Sens. Environ. 206, 287–299 (2018).ADS 

Google Scholar 
Bowman, D. M., Zhang, Y., Walsh, A. & Williams, R. Experimental comparison of four remote sensing techniques to map tropical savanna fire-scars using Landsat-TM imagery. Int. J. Wildland Fire 12, 341–348 (2003).
Google Scholar 
Hesselbarth, M. H., Sciaini, M., With, K. A., Wiegand, K. & Nowosad, J. landscapemetrics: An open-source R tool to calculate landscape metrics. Ecography 42, 1648–1657 (2019).
Google Scholar 
Koenker, R. quantreg: Quantile Regression. R package version 5.05 (R Foundation for Statistical Computing). http://CRAN.R-project.org/package=quantreg (2013).Stokes, M. A. & Smiley, T. L. Introduction to Tree-Ring Dating (University of Chicago Press, 1968).
Google Scholar 
Geoscience Australia. Surface Geology of Australia 1:1 Million Scale Dataset 2012 Edition. https://data.gov.au/data/dataset/surface-geology-of-australia-1-1-million-scale-dataset-2012-edition. More

Data UsedWe used cadastral maps from 1858 to reconstruct the LULC structure of Aksu and Kestel for the mid-nineteenth century. General Staff of the Ottoman Army produced these maps in 1:10,000 scale. These maps were one of the earliest attempts of creating cadastral maps in the Ottoman Empire. The images of historical maps scanned at 1270 dpi resolutions are provided by the Turkish Presidency State Archives of the Republic of Turkey – Department of Ottoman Archives (Map collection, HRT.h, 561–567). Individual tiles of cadastral maps are of a 1:2,000 scale. However, these maps are kept separated from their accompanying cadastral registers or documentation regarding their production process in the archives. There are no studies on the production of these cadastral maps, but few studies used them until now35,36.The LULC structures of Aksu and Kestel for the mid-twentieth century were generated using aerial photographs from June 23, 1955, with a scale of 1:30,000. These aerial photographs were captured by the US Navy Photographic Squadron VJ-62 (established on April 10, 1952, re-designated to VAP-62 on July 1956, and disestablished on October 15, 1969). The squadron conducted mapping and special photographic projects worldwide37. Lastly, the VHR satellite images of WorldView-3 (WV-3) with 0.3 m of spatial resolution, collected on September 6, 2020, were used to show the up-to-date LULC patterns of Aksu and Kestel.MethodologyFigure 2 shows the flowchart of steps followed in this study to detect the LULC changes. The workflow includes three phases: preprocessing, LULC mapping, and statistical analysis of LULC changes.Figure 2Flowchart of the processing steps for the LULC change analysis for Kestel.Full size imageData preprocessingOrthorectification is the first important step in ensuring accurate spatial positioning among the multi-temporal and multi-source images, eliminating geometric distortions, and defining all data sets on a common projection system. To align the multi-modal geospatial datasets, we first performed the orthorectification of the satellite images and then we used the orthorectified satellite images as reference for the georeferencing of the cadastral maps and aerial photographs.Satellite imagery orthorectificationWe first pan-sharpened the WV-3 images by applying the PANSHARP2 algorithm38 to fuse the panchromatic (PAN) image of 0.3 m spatial resolution with four multispectral bands (R, G, B, and near-infrared (NIR)) of 1.2 m. We then geometrically corrected the pan-sharpened (PSP) WV-3 imageries using an ALOS Global Digital Surface Model with a horizontal resolution of approximately 30 m (ALOS World 3D – 30 m), rational polynomial coefficients (RPC) file, and additional five ground control points (GCPs) for the refinement. As a geometric model, we used the RPC model with zero-order polynomial adjustment39, and orthorectified images were registered to the Universal Traverse Mercator (UTM) Zone 35 N as the reference coordinate system.Georeferencing of scanned cadastral maps and aerial photographsWe used orthorectified WV-3 imageries as a reference for the geometric correction of the historical cadastral maps and the aerial photographs. We selected the spline (triangulation) transformation, a rubber sheeting method, useful for local accuracy and requiring a minimum of 10 control points, as the transformation method to determine the correct map coordinate location for each cell in the historical maps and aerial photographs. The spline transformation provides superior accuracies for the geometric correction of the historical geospatial data40,41.The lack of topographic properties of planimetric features in the historical cadastral maps and the inherent distortions of the aerial photographs due to terrain and camera tilts causes difficulties in precise georeferencing of these data sets. It increases the number of required ground control points (GCPs) for optimal image rectification. Adequate and homogenously distributed GCPs, identified from cadastral maps and aerial photographs, can ensure precise spatial alignment among different geospatial data. The best locations for GCPs were border intersections of agricultural fields and roads. As for artificial objects, places of worship and schools, which are highly probable that have remained unchanged, are other optimal locations for GCPs to perform the accurate geometric correction. Figure 3 displays samples of GCPs selected from cadastral maps and aerial photographs. We obtained 2.11 m or better overall RMSE (Root Mean Square Error) values for the geometric correction of the historical maps and aerial photographs.Figure 3Examples of GCPs selection (red crosses in blue circles) on (a), (c) Cadastral maps and their counterparts on (b), (d) Aerial photographs.Full size imageLULC classification schemeWe defined our classification scheme by analyzing the LULC classes distinguished in all three datasets (i.e., cadastral maps, aerial photographs, and WV-3 imageries). We used the classification scheme shown in Table 1. We also present codes and names of the land cover (LC) classes derived from Corine LC nomenclature42.Table 1 Classification scheme of the study.Full size tableThe legends provided on the historical cadastral maps of Aksu and Kestel delineate 15 LULC categories, which are: (1) buildings, (2) home gardens, (3) roads, (4) arable land, (5) gardens, (6) mulberry groves, (7) chestnut groves, (8) olive groves, (9) vegetable gardens, (10) forest, (11) courtyards, (12) vineyards, (13) arable fields, (14) cemeteries, (15) watercourses. Categorizing the land cover types of cadastral maps is limited with the classes available in the map legend. The legend of cadastral maps categorizes the forested area in one class named “forest”. Therefore, it was not possible to use third-level LC sub-categories in our classification schema for forest area which is separating forested areas into three subclasses (3.1.1, 3.1.2, and 3.1.3) according to the type of tree cover. Although some of the third-level LC sub-categories could be extracted from the cadastral map legend, it was not possible to extract all third level agricultural classes from single-band monochromatic aerial photographs. Although the spatial extent of fruit trees as a permanent crop could be determined from aerial photographs, it was not possible to classify these trees into different fruit types (e.g. 2.2.1 Vineyards, 2.2.2 Fruit trees and berry plantations, 2.2.3 Olive groves). Limitation on the number of forest classes is due to the historical cadastral map content; whereas limitation on the number of agricultural classes is mainly offset by the aerial photographs which have only one spectral band and we did not have any field survey or ancillary geographical data that could help the specific identification of fruit trees.Our primary focus is to find out agricultural land abandonment, deforestation/afforestation, urbanization, and rural depopulation within the historical periods. Therefore, most of the second level LULC classes are sufficient for our purpose. LULC changes within the third class level such as the conversion of third level agriculture classes among each other were not aimed to analyze in this research. Our datasets allow us to use Level 3 CORINE classes for the artificial surfaces. These classes are useful to determine residential area implications of rural depopulation or urbanization, one of the focused transformation types for our analysis.We re-organized and categorized the LULC types used in the cadastral maps, with minimum possible manipulation (only for 2.4 and 3.2 LC classes) according to the classification scheme, as shown in Table 2.Table 2 Correspondence between Corine Land Cover and historical cadastral maps nomenclature.Full size tableLULC mappingAfter aligning all geospatial data, we used the georeferenced cadastral maps, aerial photographs, and satellite images for the LULC mapping. We set the spatial extent of the selected regions based on boundaries digitized from the cadastral maps of 1858. Then we detected historical LULC changes within these extents for all geospatial datasets covering 1900 ha and 830 ha of the Aksu and Kestel regions, respectively. Figures 4 and 5 show the selected extents from the historical maps, aerial photographs, and satellite images of the Kestel and Aksu sites, respectively.Figure 4Geospatial dataset for the Kestel study region. (a) 1858 Cadastral map, (b) 1955 aerial photo, and (c) 2020 WV-3 satellite image (finer details shown in the inset images highlighted by Blue boxes).Full size imageFigure 5Geospatial dataset for the Aksu study region. (a) 1858 Cadastral map, (b) 1955 aerial photo, and (c) 2020 WV-3 satellite image (finer details shown in the inset images highlighted by red boxes).Full size imageDigitization of cadastral maps-1858 LULC mapsWe visually interpreted and manually digitized the geographic features on the historical maps and created vector data for each class. The road features in cadastral maps lack topological properties. They also include spatial errors possibly generated in the process of surveying and map production. Therefore, we cross-checked digitized road segments by visual inspection of the road data of the aerial photographs from 1955. We then further modified road polygons to represent accurate road widths. Afterward, we categorized vectorized features of the cadastral maps into the LULC classes defined in Table 1. Finally, we created the vectorized 1858 LULC map. Figure 6 presents the vectorized 1858 cadastral maps of Aksu and Kestel.Figure 6Vectorized cadastral maps of (a) Kestel and (b) Aksu with Red and green lines showing the vector boundaries.Full size imageObject-based image analysis of aerial photographs-1955 LULC mapsAt the second stage of LULC mapping, we performed the segmentation and classification of the aerial photographs using an object-based approach for generating the 1955 LULC map. The object-based image analysis (OBIA) approach in LULC mapping provides advantages over the traditional per-pixel techniques such as higher classification accuracy, depicting more accurate LULC change, and differentiating extra LULC classes33,43,44. We used the eCognition® software (Trimble Germany GmbH, Munich) to implement an object-based image analysis (OBIA). The OBIA approach contains two phases including the segmentation and classification phases that are performed to locate meaningful objects in an image and categorize the created objects, respectively.Multiple ancillary datasets have been used to support different phases of OBIA. The Open Street Map (OSM) vector data, an open-source geospatial dataset (http://www.openstreetmap.org/), has been utilized as ancillary vector data in OBIA to improve the classification of the remotely sensed images. Sertel et al. (2018) used OSM as a thematic layer for road extraction7. Since there are several limitations in extracting the roads from aerial imagery, the OSM road network data could be useful. A majority of unpaved roads in single-band aerial photographs can easily be misclassified as homogeneous areas of arable lands. Precise detection of the roads from monoband aerial photographs without multi-spectral information is difficult. Therefore, we overlaid the OSM road network data with the aerial photographs to extract the revised aerial road vectors through visual interpretation and manual digitization.We segmented the 1955 aerial photographs with the integration of 1858 LULC map produced from cadastral maps. We implemented the multi-resolution segmentation algorithm. In this segmentation method, a parameter called scale determines the size of resulting objects, and the shape and compactness parameters determine the boundaries of objects. The segmentation process of the aerial photographs was performed at multiple stages with various scale, shape, and compactness parameter values. At the initial stage, we segmented the regions according to the 1858 LULC map and we utilized large-scale parameters. The scale parameter was set to 100 and the shape parameter and the compactness were set as 0.7 and 0.3, respectively. At this stage, we focused on interpreting the objects that have not changed between 1858 and 1955. We classified the segments using the thematic layer attribute (LULC classes defined by the cadastral maps) with the highest coverage. Image objects in which the land surface has changed during 1858–1955 period were detected by visual interpretation and unclassified for further segmentation. We followed this approach to reduce the manual effort. We defined unchanged objects between 1858 and 1955 and assigned the same classes of 1858 LULC map to the objects in 1955 aerial photographs. We then segmented the remaining segments, the last time into smaller objects with the scale parameter set as 25, the shape parameter set as 0.2, and the compactness set as 0.8.We classified the remaining unclassified objects through the development of rulesets. An object can be described by several possible features as explanatory variables which are provided by eCognition. In the classification ruleset, different features and parameters can be defined to describe and extract object classes of interest and thresholds for each feature can be defined by the trial-and-error method. We tested sets of variables for the classification of the monoband aerial photographs. Object features such as the mean value of the monoband, texture after Haralick, distance to neighbor objects, shape features (e.g., rectangular fit and asymmetry), and extent features (e.g., area and length/width) were the most useful alternatives. The classification process of the parcels of the aerial photographs with LULC change started with the classification of roads constructed between 1858 and 1955 by utilizing the aerial road map. The watercourse class was the most difficult to classify since shrubs or trees mostly covered the watercourses. These areas were misclassified as forest or agricultural land. Therefore, experts in historical map reading with local geographical information performed the detection and classification of the water course class and interpreted by the cadastral map (1858) and the google map (2020). After roads and watercourses, we classified forest and agricultural lands using the optimal thresholds for the brightness feature. We calculated the thresholds using the single band of the aerial photograph combined with the area and rectangular fit features. The heterogeneous agricultural areas class principally occupied by agriculture with significant areas of natural grass and trees within the same object are separated from the arable lands using the standard deviation of the digital number (DN) values of the aerial photographs. The texture feature helped classify the permanent crops. The brightness, shape, asymmetry, and distance to road class features were the best-performing ones for classifying the remaining artificial surfaces. The manual interpretation was performed for the classification of sub-classes of artificial surface class, including the continuous/discontinuous urban fabric, industrial, commercial, and transport units, mine, dump and construction sites, and artificial, non-agricultural vegetated areas. Since these land use classes contain one or more land cover and land use categories (e.g., artificial non-agriculture land or industrial or commercial units), finding the optimal threshold and exact feature for distinguishing the subclasses of artificial surfaces is difficult. Especially in the case of using the single-band aerial photographs, manual interpretation was required.Object-based image analysis of satellite images-2020 LULC mapsWe segmented WV-3 satellite images using multi-resolution segmentation algorithm and ancillary geographic data. Similar to the aerial road map, the road network of the study region in 2020, named, WV-3 road map, was extracted by overlaying the OSM road data with the WV-3 satellite image. In the segmentation process of the WV-3 image, we used the vector boundaries of the classified aerial photograph (the 1955 LULC map) and the WV-3 road map as ancillary thematic layers. We opted for the same segmentation and classification approach used for the aerial photographs for the WV-3 image.Firstly, we segmented the satellite image into spectrally homogeneous objects using vector data of the 1955 LULC map by applying large-scale parameters. We implemented scale parameter values of 300, 200, 100, and 50 to find the optimal scale to classify objects that have not changed between 1955 and 2020. The best multi-resolution segmentation configuration was the scale of 100 and the shape and compactness parameters of 0.3 and 0.7, respectively. We classified the segments using the thematic layer attribute (LULC classes defined by the aerial maps) with the highest coverage. Segments with LULC change, e.g. the image objects in which the land surface has changed during 1955–2020 period were detected by visual interpretation and unclassified for further segmentation. As a result, we excluded the objects which were remained unchanged during 1955–2020 by assigning the prepared labels which were allocated in the previous step during the classification of 1955 aerial photographs. We then segmented the remaining objects into smaller objects to identify the changed areas in detail. At this step, the scale, shape, and compactness parameters were set as 25, 0.2, and 0.8, respectively.Except for the additional sets of variables utilized to classify the WV-3 images, we applied the rule-set developed for the classification of the aerial photograph for the classification of the remaining objects of 2020 satellite images. The additional sets of variables include the mean of G, B, R, and NIR and two spectral indices, the Normalized Difference Water Index (NDWI), and the Normalized Difference Vegetation Index (NDVI). NDVI was calculated as the normalized difference of reflectance values in the red and NIR bands; whereas , NDWI was determined as the normalized difference of reflectance values of the green and NIR bands. Through the logical conditions, objects having specified values of NDVI and NDWI can be assigned to vegetation and water classes, respectively. The use of NDVI facilitated the delineation of terrains covered by vegetation and the NDWI improved the extraction of water bodies due to its ability to separate water and non-water objects. We separated different sub-classes of agricultural areas and forests by using optimal thresholds for NDVI which were defined by a trial and error method. Also we utilized assigning the optimal threshold to NDWI to separate water bodies from other land covers. In addition, the mean blue band layer was useful in classifying the artificial surfaces. We assessed the accuracy of each classification using error matrices (overall, user’s and producer’s accuracies, and Kappa statistics)45,46.Estimating LULC changes and LULC conversionsAfter the production of LULC maps of Aksu and Kestel for 1858, 1955, and 2020, the vector data of the LULC maps were used to quantify the LULC conversions for two different periods which are 1858–1955 and 1955–2020. To compare the LULC maps of study areas between two different dates of each study period, we provided detailed “from-to” LULC change information by calculating the LULC change transition matrix computed using overlay functions in ArcGIS.We overlaid LULC maps of 1858 and 1955 and intersected the vector boundaries of the 1858 and 1955 LULC maps to determine the conversion types of LULC classes (from which class to which class). Similarly, to quantify the LULC changes between 1955 and 2020, we overlaid the 1955 and 2020 LULC maps. Then we created transition matrices and performed statistical analysis utilizing the matrices. Finally, we discussed the main LULC change types and the driving factors of the changes in the selected study areas. More

Drainage area and branching ratio: a matter of scaleGeomorphological and ecological viewpoints on river networks generally differ owing to discordant definitions of the fundamental unit (the node) used to analyze them. From a geomorphological perspective, the determination of a river network entails the definition of an observational scale. Real river networks can be extracted from digital elevation models (DEMs) via algorithms for flow direction determination such as D8 (i.e., each pixel drains towards the lowest of its 8 nearest neighbors53). After the outlet location has been specified (and hence the upstream area A spanned by the river network), the first observational scale required is thus the pixel length l of the DEM, which defines the extent of a network node. A second scale is then needed to distinguish the portion of the drainage network effectively belonging to the channel network. The simplest but still widely used method53 defines channels as those pixels whose drainage area exceeds a threshold value AT. Hydrologically based criteria to determine the appropriate value for AT exist54; however, for the sake of simplicity, we here consider AT as a free parameter.BBTs and RBNs are random constructs, and as such they do not satisfy the optimality criterion of minimizing total energy expenditure, which is the fundamental physical process shaping fluvial landscapes. Furthermore, neither of these networks is a spanning tree, which is a key attribute of real fluvial landforms10: in fact, in both BBTs and RBNs, the extent of the drained domain is not defined. As a result, the drainage area at an arbitrary network node cannot in principle be attributed, unless by using the number of upstream nodes as a proxy. This has practical implications from an ecological viewpoint because drainage area is the master variable controlling several attributes of a river, such as width, depth, discharge, or slope3,55, which in turn impact habitat characteristics and the ecology of organisms therein56.In BBTs and RBNs, branching probability p has been defined35,38,45,46,47 as the probability that a network node is branching, i.e. connected to two upstream nodes. As such, the branching probability of a realized river network (be it a real river or a synthetic construct) could be evaluated as the ratio between the number of links NL constituting a network and the total number of network nodes N; if a unit distance between two adjacent nodes is assumed, the denominator equals the total network length. We note that the former definition of branching probability only holds in the context of the generation of a synthetic random network; it is in fact improper to refer to a “probability” when analyzing the properties of a realized river network. We clarify this aspect by introducing the concept of branching ratio pr for the latter definition (pr = NL/N). Moreover, in the case of BBTs, p and pr do not coincide (see Methods). Importantly, p and pr have no parallel in the literature on fluvial forms, nor do they refer to any of the well-studied measures of rivers’ fractal character.The choice of different observational scales for the same drainage network results in different values of NL and N, and hence of pr. Remarkably, the very same drainage network can result in river networks that virtually assume any value of pr (ranging from 0 to 1) and N (up to the upper bound A) depending on the choice of AT and A (the latter corresponding to a given l value when measured in the number of pixels; Fig. 1d–i); networks with low AT/A ratios result in high N (Fig. 2a), while networks with low AT result in high pr (Fig. 2b). Furthermore, pr does not identify the inherent (i.e., scale-independent) branching character of a given river network in relation to other river networks. In fact, by extracting different river networks at various scales (i.e., various AT values) and assessing the rivers’ rank in terms of pr, one observes that rivers that look more “branching” (i.e., have higher pr) than others for a given AT value can become less “branching” for a different AT value (Fig. 3). We therefore conclude that branching probability is a non-descriptive property of a river network, which by no means describes its inherent branching character, and depends on the observational scale.Fig. 2: Variation of N and pr as a function of observational scales for OCNs and real river networks.a Expected value of number of network nodes N as a function of threshold area AT and total drained area A (from Eq. (1)); the white dots indicate the values of AT and A used to generate the OCNs used in this analysis. b Expected value of branching ratio pr as a function of AT and A (from Eq. (1)); symbols as in a.Full size imageFig. 3: Values of branching ratio as a function of AT for the 50 real river networks analyzed in this study.a Natural values of pr in logarithmic scale. b z-normalized branching ratios (i.e., for each AT value, values of pr are normalized so that they have null mean and unit standard deviation), which better shows how rivers rank differently in terms of pr for different observation scales (i.e., AT). Lines connect dots relative to the same river. For visual purposes, rivers that rank first, second, second-to-last or last in at least one of the AT groups are displayed in colors; the other rivers are displayed in grey.Full size imageScaling is also crucial when looking at river networks from an ecological perspective. In this case, the relevant scale determining the dimension l of a node is the extent of habitat within which individuals (due to e.g. physical constrains) can be assigned to a single population57,58; the riverine connectivity and ensuing dispersal among these populations give rise to a metapopulation at the river network level. The specific spatial scale largely depends on the targeted species (e.g. being larger for fish than for aquatic insects), and it is conceivably much larger than (or, at least, it has no reason to be equal to) the pixel size of the DEM on which the river network is extracted. Since the evaluation of pr depends on the number of nodes N, which, in turn, is defined based on the scale length l, the resulting pr of a river network under this perspective would depend on the characteristics of the target taxa, which is inconsistent with the alleged role of pr as a scale-invariant property of river networks.Note also that using the ecological definition of l (i.e., spatial range of a local population) to discretize a real river network into N nodes, and from there calculate the branching ratio pr = NL/N, is problematic. Indeed, this would imply an elongation of all links shorter than l (which constitute a non-negligible fraction of the total links, under the assumption of exponential distribution of link lengths51), hence preventing a correct estimation of the connectivity patterns (i.e., distances between nodes) and the resulting ecological metrics of the river network (see section Ecological implications).From an ecological perspective, it could be reasonable to consider AT as a parameter expressing how a particular taxon perceives the suitable landscape, rather than a value to be determined from geomorphological arguments: for instance, large fishes inhabit wide and deep river reaches, and do not access small headwaters56. In this case, imposing a large AT would result in a coarser, less branching network constituted by few main channels (Fig. 1f, i), which could mimic the potentially available habitat for such species. Conversely, aquatic insects inhabit also small headwaters17,59, therefore their perceived landscape would resemble the finely resolved networks of Fig. 1d, g, characterized by low AT and higher (apparent) pr.Topology and scaling of river networks and random analoguesTo verify the topological (i.e., Horton’s laws on bifurcation and length ratios) and scaling (i.e., probability distribution of drainage areas) relationships of the different network types, we extracted from DEMs 50 real river networks encompassing a wide range of drainage areas (Fig. 4), and we generated 50 OCNs, 50 RBNs and 50 BBTs of comparable size (see Methods).Fig. 4: Location of real river basins used in the analysis.River basins are shown in dark grey; countries in light grey. Rivers’ numbering is sorted in ascending order according to drainage area values.Full size imageTypical values3,7,60 for the bifurcation ratio RB lie between 3 and 5, while length ratios (RL) range between 1.5 and 3.5. As expected, we observed that the real rivers and OCNs used in our analysis have RB and RL values within the aforementioned ranges (Fig. 5a, b). The same is true for RBNs, while the RB and RL values found for BBTs are lower than the typical ranges. This finding holds regardless of the scale (subsumed by AT) at which real river networks and OCNs are extracted (Supplementary Figs. 1 and 2). Remarkably, BBTs fail to satisfy Horton’s laws despite the statistical inevitability of such laws for any network argued by ref. 61. To this regard, we note that the networks analyzed by ref. 61 did not include constructs where all paths from the source nodes to the outlet have the same length, which is the defining feature of BBTs (Fig. 1a).Fig. 5: Comparison of topological and scaling properties of the different networks.a Scaling of number of network links Nω as a function of stream order ω for the various network types (rivers and OCNs obtained with AT = 20 pixels; RBNs and BBTs derived accordingly – see Methods). b Mean link length Lω (in units of l) as a function of ω. Networks used are as in panel a. c Scaling of drainage areas: probability P[A ≥ a] to randomly sample a node with drainage area A ≥ a as a function of a. The displayed trend lines are fitted on the ensemble values for the 50 network replicates, by excluding nodes with drainage area larger than 2000 pixels (cutoff value marked with a black solid line). The scaling coefficients β reported correspond to the slopes of the fitted trend lines. Extended details on all panels are provided in the Supplementary Methods.Full size imageWhile the power-law scaling of areas in OCNs (Fig. 5c) has an exponent β ≈ 0.45 that closely resembles the one found for the real rivers (β ≈ 0.46) and within the typically observed range8,10 β = 0.43 ± 0.02, drainage areas of RBNs scale as a power law with an exponent β ≈ 0.51, which departs from the observed range. Conversely, BBTs do not show any power-law scaling of areas. Scaling exponents of drainage areas fitted separately for each real river network yielded values in the range 0.36÷0.57 (Supplementary Table 1). In particular, we observed that these values tend to the expected range β = 0.43 ± 0.02 for increasing values of A, expressed in number of pixels (Supplementary Fig. 3), hence implying that highly resolved catchments are required in order to properly estimate β. Interestingly, the observed values of Horton ratios and scaling exponent β for RBNs are compatible with the values RB = 4, RL = 2, β = 0.5 predicted for Shreve’s random topology model3,60,62, which is actually equivalent to a RBN with infinite links.Ecological implicationsWe compared the different network types via two metrics that express the ecological value of a landscape for a metapopulation: the coefficient of variation of a metapopulation CVM and the metapopulation capacity λM. The coefficient of variation of a metapopulation63 is a measure of metapopulation stability (a metapopulation being more stable the lower CVM is), while the metapopulation capacity42,64 expresses the potential for a metapopulation to persist in the long run (persistence being more likely the higher λM is). Both measures are among the most universal metrics describing dynamics of spatially fragmented populations24,40. In order to assess the impact of the two landscape features mostly affecting metapopulation dynamics, i.e. spatial connectivity and spatial distribution of habitat patches, we calculated these metrics for the four network types under two different scenarios: uniform (CVM,U, λM,U) and non-uniform (CVM,H, λM,H) spatial distribution of habitat patch sizes. In the first scenario, CVM,U and λM,U assess stability and persistence (respectively) of a metapopulation solely based on pairwise distances between network nodes; in the second scenario, CVM,H and λM,H depend on the interplay between pairwise distances and spatially heterogeneous habitat availability (namely, downstream nodes being larger than upstream ones).We found that the values of CVM (be it derived with uniform (CVM,U) or nonuniform (CVM,H) distributions of patch sizes) obtained for OCNs match strikingly well those of real rivers (Fig. 6). These CVM values are consistently lower than those found for RBNs, while values of CVM for BBTs are even higher. Notably, this result holds for different values of AT (and hence different pr values) at which real rivers and OCNs are extracted (Fig. 6a–c; g–i), and for values of mean dispersal distance α (see Methods) spanning multiple orders of magnitude (Supplementary Figs. 4–7).Fig. 6: Comparison of values of metapopulation metrics across river network types and observational scales (AT).a–c CVM,U. d–f λM,U. g–i CVM,H. j–l λM,H. Boxplot elements are as follows: center line, median; notches, (pm 1.58cdot {{{{{{{rm{IQR}}}}}}}}/sqrt{50}), where IQR is the interquartile range; box limits, upper and lower quartiles; whiskers, extending up to the most extreme data points that are within ±1.5 ⋅ IQR; circles, outliers. Metapopulation metric values were obtained by setting α = 100 l. Note that in Eq. (1), given A = 40, 000, AT = 20 results in E[N] ≈ 4574, E[pr] ≈ 0.228; AT = 100 yields E[N] ≈ 2231, E[pr] ≈ 0.098; AT = 500 results in E[N] ≈ 1088, E[pr] ≈ 0.042.Full size imageFor a constant α value, the CVM of real rivers, OCNs and RBNs decreases as the resolution at which the network is extracted increases (i.e., AT decreases; see Fig. 6 and Supplementary Figs. 4–7). This is expected63, since a decrease in AT corresponds to an increase in N (Fig. 2a), leading to a decrease in CVM. Indeed, a larger ecosystem, constituted of more patches, has the potential to include a larger (and more diverse) number of subpopulations, which increases stability at a metapopulation level through statistical averaging–a phenomenon widely known as the portfolio effect65. We also found that BBT networks do not generally follow the above-described pattern of decreasing CVM with increasing N; rather, the CVM of BBTs increases with N when the mean dispersal distance α is set to intermediate to high values (Fig. 6 and Supplementary Figs. 5–7), and only when α is very low (e.g. α = 10 l as in Supplementary Fig. 4) and a uniform patch-size distribution is assumed does CVM,U follow the expected decreasing trend with increasing N.However, we need to warn against the conclusion that river networks with higher values of pr (and hence lower AT, see Fig. 2b) are inherently associated with higher metapopulation stability. Indeed, our result was obtained by changing the scale at which we observed the same river networks, and not by increasing the river networks’ size. If the number of network nodes (and, consequently, the branching ratio pr) is determined by the scale at which the landscape is observed, one cannot directly assume that any of such nodes is a node (or patch) in the ecological sense, i.e. the geographical span of a local population: the extent of such patches should be determined based on the mobility characteristics of the focus species, and should be independent of the scale at which the river network is observed. In contrast, we note that, if different river networks spanning different catchment areas (say, in km2) are compared, all of them extracted from the same DEM (same l and same AT in km2), then the larger river network will appear more branching (i.e., have larger pr). Indeed, by selecting catchments with larger A (in km2) for fixed l and AT (in km2), one moves towards the top-left corner of Fig. 2a, b (i.e., perpendicular to the level curves AT/A). The apparent higher “branchiness” of the river network with larger A will result in lower values of CVM; however, the higher metapopulation stability of the larger network will not be due to its (alleged) inherent more branching character, but only dictated by its larger habitat availability.We observed that metapopulation capacity λM values of OCNs (be it evaluated under uniform (λM,U) or non-uniform (λM,H) patch-size distribution assumption) are the closest to those of real rivers, while RBNs (and even more so BBTs) generally overestimate λM with respect to real rivers and OCNs (Fig. 6d–f; j–l). This result holds irrespective of the choice of AT and for intermediate to high values of α (Supplementary Figs. 5–7). When the mean dispersal distance is instead set to very low values (α = 10 l – Supplementary Fig. 4) and the river network is extracted at a high resolution (i.e., low AT), the metapopulation capacity of OCNs under assumption of uniform patch-size distribution (λM,U) is underestimated with respect to that of real rivers. A likely explanation for this apparent mismatch is that, for low values of AT, the number of nodes N tends to be somewhat higher for the extracted river networks used in this analysis than for OCNs (Supplementary Fig. 8), and the effect of the different dimensionality of real rivers and OCNs in the metapopulation capacity estimation tends to be more evident as the mean dispersal distance decreases. Interestingly, such mismatch is absent when a non-uniform patch size distribution is assumed, as λM,H values for OCNs match those for real rivers regardless of the mean dispersal distance value and the river network resolution (Fig. 6; Supplementary Figs. 4–7).The OCN construct encapsulates both random and deterministic processes, the former related to the stochastic nature of the OCN generation algorithm, and the latter pertaining to the minimization of total energy expenditure that characterizes OCN configurations. As such, OCNs reproduce the aggregation patterns of real river networks. From an ecological viewpoint, this implies that both pairwise distances between nodes and the distribution of patch sizes (expressed as a function of drainage areas, or of a proxy thereof such as the number of nodes upstream) are much closer to those of real networks than is the case for fully random synthetic networks as BBTs and RBNs. In particular, BBTs and (to a lesser extent) RBNs tend to underestimate pairwise distances with respect to real rivers and OCNs, as documented by a comparison of mean pairwise distances across network types (Supplementary Fig. 9a–c). Our analysis shows that the connectivity structure of these random networks (subsumed by the matrix of pairwise distances) is too compact with respect to that of real rivers, which leads to an overestimation of the role of dispersal in increasing the ability of a metapopulation to persist in the long run, but also an increased likelihood of synchrony among the different local populations, which results in higher instability.Comparison of patch size distributions among the network types expressed in terms of CVM,0 (i.e., the portion of CVM,H that uniquely depends on the distribution of patch sizes and not on pairwise distances) shows that, while for coarsely resolved networks (AT = 500) no clear differences in CVM,0 emerged, for highly resolved networks (AT = 20) BBTs heavily underestimate the CVM,0 of real rivers and OCNs, while RBNs slightly overestimate it (Supplementary Fig. 9d–f). As a result of the interplay of differences in distance matrices and patch size distributions, BBTs and (to a lesser extent) RBNs generally tend to overestimate the coefficient of variation of a metapopulation and the metapopulation capacity of real rivers and OCNs in both scenarios of uniform and non-uniform patch size distribution. The only exception to this trend occurs for the metapopulation capacity λM,H of very large BBTs (corresponding to AT = 20) in the case of very high dispersal distances (α = 1000 l – Supplementary Fig. 7): here, the patch-size effect (i.e., underestimation of CVM,0) predominates over the distance effect (i.e., overestimation of mean dij), resulting in an underestimation of λM,H with respect to real rivers and OCNs.Our results were derived under a number of simplifying assumptions. In particular, we acknowledge that, while the distance matrix of a landscape and the distribution of patch sizes have in general important implications for metapopulation dynamics, other factors not considered here, such as Euclidean between-patch distance48, fat-tailed dispersal kernel66 and density-dependent dispersal67 could also play a relevant role in this respect. However, it needs to be noted that, especially with regards to the assessment of the Moran effect in metapopulation synchrony (i.e., increased synchrony in local fluvial populations that are geographically close but not flow-connected48), the use of OCNs allows integration of Euclidean distances in a metapopulation model, while this is not possible for RBNs and BBTs, where Euclidean distances are not defined. Moreover, if a larger degree of realism is required for a specific ecological modelling study, such as heterogeneity in abiotic factors (e.g. water temperature or flow rates), the use of OCNs as model landscapes allows a direct integration of these variables, as they can conveniently be expressed as functions of drainage area3,55. In contrast, this is not possible for RBNs or BBTs, because only OCNs verify the scaling of areas (Fig. 5c), while RBNs and BBTs lack a proper definition of drainage areas.Our comparison of synthetic and real river networks showed that riverine metapopulations are more stable and less invasible than what would be predicted by random network analogues. Conversely, the use of OCNs as model landscapes allows capturing not only the scaling features of real rivers, but also drawing ecological conclusions that are in line with those that could be observed in real river networks. We thus support the use of OCNs as analogues of real river networks in theoretical and applied ecological modelling studies. While we found that BBTs are highly inaccurate in reproducing ecological metrics of real river networks and should be therefore discarded altogether in future modelling applications, RBNs show a certain degree of similarity with OCNs and real river networks in this respect; moreover, RBNs (as is the case for any random tree61) satisfy Horton’s laws on bifurcation and length ratios. A relevant advantage of RBNs over OCNs is that their generation algorithm is at least one order of magnitude faster49. Therefore, we acknowledge that RBNs could be considered as a suitable surrogate for real river networks as null models in cases where a large number of network replicates is required. To this end, we encourage researchers exploiting synthetic river networks (whether they be OCNs or RBNs) to always clarify the observational scales (that is, total area drained, size of a node, area drained by a headwater) subsumed by the synthetic network and which give rise to a certain complexity measure (i.e., branching ratio). Only in such a way could the predictions from these studies be compared with real river networks.In conclusion, our results advocate a tighter integration between physical (geomorphology, hydrology) and biological (ecology) disciplines in the study of freshwater ecosystems, and particularly in the perspective of a mechanistic understanding of drivers of persistence and loss of biodiversity. More

Harvey, P. H. & Pagel, M. D. The comparative method in evolutionary biology. Vol. 239 (Oxford University Press, 1991).Oyston, J. W., Hughes, M., Wagner, P. J., Gerber, S. & Wills, M. A. What limits the morphological disparity of clades? Interface Focus 5, 0042 (2015).Article 

Google Scholar 
Jetz, W., Thomas, G. H., Joy, J. B., Hartmann, K. & Mooers, A. O. The global diversity of birds in space and time. Nature 491, 444–448 (2012).CAS 
PubMed 
Article 

Google Scholar 
Webb, C. O. Exploring the phylogenetic structure of ecological communities: an example for rain forest trees. Am. Naturalist 156, 145–155 (2000).Article 

Google Scholar 
Purvis, A., Gittleman, J. L. & Brooks, T. Phylogeny and conservation. (Cambridge University Press, 2005).Page, R. D. M. Parallel phylogenies: reconstructing the history of host-parasite assemblages. Cladistics 10, 155–173 (1994).Article 

Google Scholar 
Weaver, S. C. & Vasilakis, N. Molecular evolution of dengue viruses: contributions of phylogenetics to understanding the history and epidemiology of the preeminent arboviral disease. Infect., Genet. Evolution 9, 523–540 (2009).CAS 
Article 

Google Scholar 
Tassy, P. Trees before and after Darwin. J. Zool. Syst. Evolut. Res. 49, 89–101 (2011).Article 

Google Scholar 
Heather, J. M. & Chain, B. The sequence of sequencers: The history of sequencing DNA. Genomics 107, 1–8 (2016).CAS 
PubMed 
Article 

Google Scholar 
Pyron, R. A. Post-molecular systematics and the future of phylogenetics. Trends Ecol. Evolution 30, 384–389 (2015).Article 

Google Scholar 
Sansom, R. S. & Wills, M. A. Differences between hard and soft phylogenetic data. Proc. R. Soc. B: Biol. Sci. 284, 20172150 (2017).Article 

Google Scholar 
Scotland, R. W., Olmstead, R. G. & Bennett, J. R. Phylogeny reconstruction: the role of morphology. Syst. Biol. 52, 539–548 (2003).PubMed 
Article 

Google Scholar 
Regier, J. C. et al. Arthropod relationships revealed by phylogenomic analysis of nuclear protein-coding sequences. Nature 463, 1079–1083 (2010).CAS 
PubMed 
Article 

Google Scholar 
Callender-Crowe, L. M. & Sansom, R. S. Osteological characters of birds and reptiles are more congruent with molecular phylogenies than soft characters are. Zool. J. Linn. Soc. 194, 1–13 (2022).Article 

Google Scholar 
Wahlberg, N. et al. Synergistic effects of combining morphological and molecular data in resolving the phylogeny of butterflies and skippers. Proc. R. Soc. B: Biol. Sci. 272, 1577–1586 (2005).CAS 
Article 

Google Scholar 
He, L. et al. A molecular phylogeny of selligueoid ferns (Polypodiaceae): Implications for a natural delimitation despite homoplasy and rapid radiation. Taxon 67, 237–249 (2018).Article 

Google Scholar 
Fernández, R., Edgecombe, G. D. & Giribet, G. Phylogenomics illuminates the backbone of the Myriapoda Tree of Life and reconciles morphological and molecular phylogenies. Sci. Rep. 8, 1–7 (2018).
Google Scholar 
Eme, L., Spang, A., Lombard, J., Stairs, C. W. & Ettema, T. J. G. Archaea and the origin of eukaryotes. Nat. Rev. Microbiol. 15, 711–723 (2017).CAS 
PubMed 
Article 

Google Scholar 
Asher, R. J., Bennett, N. & Lehmann, T. The new framework for understanding placental mammal evolution. BioEssays 31, 853–864 (2009).CAS 
PubMed 
Article 

Google Scholar 
Shoshani, J. & McKenna, M. C. Higher taxonomic relationships among extant mammals based on morphology, with selected comparisons of results from molecular data. Mol. Phylogenetics Evolution 9, 572–584 (1998).CAS 
Article 

Google Scholar 
Beck, R. M. D. & Baillie, C. Improvements in the fossil record may largely resolve current conflicts between morphological and molecular estimates of mammal phylogeny. Proc. R. Soc. B: Biol. Sci. 285, 20181632 (2018).Article 

Google Scholar 
Zou, Z. T. & Zhang, J. Z. Morphological and molecular convergences in mammalian phylogenetics. Nat. Commun. 7, 1–9 (2016).
Google Scholar 
Hillis, D. M. Molecular versus morphological approaches to systematics. Annu. Rev. Ecol. Syst. 18, 23–42 (1987).Article 

Google Scholar 
Thompson, N. Alfred Russell Wallace Contributions to the theory of Natural Selection, 1870, and Charles Darwin and Alfred Wallace, ‘On the Tendency of Species to form Varieties’ (Papers presented to the Linnean Society 30th June 1858). (Routledge, 2004).Croizat, L. Panbiogeography; or an introductory synthesis of zoogeography, phytogeography, and geology, with notes on evolution, systematics, ecology, anthropology, etc., Vol. 1, 2a & 2b (Published by the author, Caracas., 1958).Means, J. C. & Marek, P. E. Is geography an accurate predictor of evolutionary history in the millipede family Xystodesmidae? PeerJ 5, e3854 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
Wills, M. A., Barrett, P. M. & Heathcote, J. F. The modified gap excess ratio (GER*) and the stratigraphic congruence of dinosaur phylogenies. Syst. Biol. 57, 891–904 (2008).PubMed 
Article 

Google Scholar 
Fisher, D. C. Stratocladistics: integrating temporal data and character data in phylogenetic inference. Annu. Rev. Ecol., Evolution Syst. 39, 365–385 (2008).Article 

Google Scholar 
Lazarus, D. B. & Prothero, D. R. The role of stratigraphic and morphologic data in phylogeny. J. Paleontol. 58, 163–172 (1984).
Google Scholar 
Camerini, J. R. Evolution, biogeography, and maps: an early history of Wallace’s Line. Isis 84, 700–727 (1993).CAS 
PubMed 
Article 

Google Scholar 
Upchurch, P., Hunn, C. A. & Norman, D. B. An analysis of dinosaurian biogeography: evidence for the existence of vicariance and dispersal patterns caused by geological events. Proc. R. Soc. B: Biol. Sci. 269, 613–621 (2002).Article 

Google Scholar 
Ferreira, G. S., Bronzati, M., Langer, M. C. & Sterli, J. Phylogeny, biogeography and diversification patterns of side-necked turtles (Testudines: Pleurodira). R. Soc. Open Sci. 5, 171773 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Ronquist, F. & Sanmartín, I. Phylogenetic methods in biogeography. Annu. Rev. Ecol., Evolution, Syst. 42, 441–464 (2011).Article 

Google Scholar 
IUCN. The IUCN Red List of Threatened Species. Version 2019-2., https://www.iucnredlist.org (2019).GBIF.org. GBIF Home Page, https://www.gbif.org/ (2019).Uetz, P., Freed, P., Aguilar, R. & Hošek, J. The reptile database., http://www.reptiledatabase.org (2019).Archie, J. W. Homoplasy excess ratios: new indices for measuring levels of homoplasy in phylogenetic systematics and a critique of the consistency index. Syst. Zool. 38, 253–269 (1989).Article 

Google Scholar 
Wilkinson, M. On phylogenetic relationships within Dendrotriton (Amphibia: Caudata: Plethodontidae) is there sufficient evidence? Herpetological J. 7, 55–65 (1997).
Google Scholar 
O’Connor, A. & Wills, M. A. Measuring stratigraphic congruence across trees, higher taxa, and time. Syst. Biol. 65, 792–811 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
Colless, D. H. Review of phylogenetics: the theory and practice of phylogenetic systematics. Syst. Zool. 31, 100–104 (1982).Article 

Google Scholar 
Lartillot, N. & Philippe, H. Improvement of molecular phylogenetic inference and the phylogeny of Bilateria. Philos. Trans. R. Soc. B: Biol. Sci. 363, 1463–1472 (2008).Article 

Google Scholar 
Sansom, R. S., Choate, P. G., Keating, J. N. & Randle, E. Parsimony, not Bayesian analysis, recovers more stratigraphically congruent phylogenetic trees. Biol. Lett. 14, 20180263 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Rosa, B. B., Melo, G. A. & Barbeitos, M. S. Homoplasy-based partitioning outperforms alternatives in Bayesian analysis of discrete morphological data. Syst. Biol. 68, 657–671 (2019).PubMed 
Article 

Google Scholar 
Lucena, D. A. & Almeida, E. A. Morphology and Bayesian tip-dating recover deep Cretaceous-age divergences among major chrysidid lineages (Hymenoptera: Chrysididae). Zool. J. Linn. Soc. 194, 36–79 (2022).Article 

Google Scholar 
O’Reilly, J. E. et al. Bayesian methods outperform parsimony but at the expense of precision in the estimation of phylogeny from discrete morphological data. Biol. Lett. 12, 20160081 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
Smith, M. R. Bayesian and parsimony approaches reconstruct informative trees from simulated morphological datasets. Biol. Lett. 15, 20180632 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
Wiens, J. The role of morphological data in phylogeny reconstruction. Syst. Biol. 53, 653–661 (2004).PubMed 
Article 

Google Scholar 
O’Leary, M. A. & Kaufman, S. G. MorphoBank 3.0: Web application for morphological phylogenetics and taxonomy., http://www.morphobank.org (2012).de Queiroz, A. & Gatesy, J. The supermatrix approach to systematics. Trends Ecol. Evolution 22, 34–41 (2007).Article 

Google Scholar 
Wilkinson, M. A comparison of two methods of character construction. Cladistics 11, 297–308 (1995).Article 

Google Scholar 
Brazeau, M. D. Problematic character coding methods in morphology and their effects. Biol. J. Linn. Soc. 104, 489–498 (2011).Article 

Google Scholar 
Drummond, A. J., Ho, S. Y. W., Phillips, M. J. & Rambaut, A. Relaxed phylogenetics and dating with confidence. PLoS Biol. 4, e88 (2006).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
O’Reilly, J. E., Puttick, M. N., Pisani, D. & Donoghue, P. C. Probabilistic methods surpass parsimony when assessing clade support in phylogenetic analyses of discrete morphological data. Palaeontology 61, 105–118 (2018).PubMed 
Article 

Google Scholar 
Keating, J. N., Sansom, R. S., Sutton, M. D., Knight, C. G. & Garwood, R. J. Morphological phylogenetics evaluated using novel evolutionary simulations. Syst. Biol. 69, 897–912 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Makarenkov, V. et al. Weighted bootstrapping: a correction method for assessing the robustness of phylogenetic trees. BMC Evolut. Biol. 10, 1–16 (2010).Article 
CAS 

Google Scholar 
Stayton, C. T. The definition, recognition, and interpretation of convergent evolution, and two new measures for quantifying and assessing the significance of convergence. Evolution 69, 2140–2153 (2015).PubMed 
Article 

Google Scholar 
Sattler, R. Homology – a continuing challenge. Syst. Bot. 9, 382–394 (1984).Article 

Google Scholar 
Jenner, R. A. & Schram, F. R. The grand game of metazoan phylogeny: rules and strategies. Biol. Rev. 74, 121–142 (1999).Article 

Google Scholar 
Pisani, D. & Wilkinson, M. Matrix representation with parsimony, taxonomic congruence, and total evidence. Syst. Biol. 51, 151–155 (2002).PubMed 
Article 

Google Scholar 
Arcila, D. et al. Testing the utility of alternative metrics of branch support to address the ancient evolutionary radiation of tunas, stromateoids, and allies (Teleostei: Pelagiaria). Syst. Biol. 70, 1123–1144 (2021).PubMed 
Article 

Google Scholar 
Felsenstein, J. Phylogenies and the comparative method. Am. Naturalist 125, 1–15 (1985).Article 

Google Scholar 
Bremer, K. Branch support and tree stability. Cladistics 10, 295–304 (1994).Article 

Google Scholar 
Johnson, W. E. et al. The late Miocene radiation of modern Felidae: a genetic assessment. Science 311, 73–77 (2006).CAS 
PubMed 
Article 

Google Scholar 
Van der Made, J. Biogeography and climatic change as a context to human dispersal out of Africa and within Eurasia. Quat. Sci. Rev. 30, 1353–1367 (2011).Article 

Google Scholar 
May, F., Rosenbaum, B., Schurr, F. M. & Chase, J. M. The geometry of habitat fragmentation: Effects of species distribution patterns on extinction risk due to habitat conversion. Ecol. Evolution 9, 2775–2790 (2019).Article 

Google Scholar 
Swofford, D. L. et al. Bias in phylogenetic estimation and its relevance to the choice between parsimony and likelihood methods. Syst. Biol. 50, 525–539 (2001).CAS 
PubMed 
Article 

Google Scholar 
Jaeger, J. J. & Martin, M. African marsupials – vicariance or dispersion? Nature 312, 379–379 (1984).Article 

Google Scholar 
Smith, B. T. et al. The drivers of tropical speciation. Nature 515, 406–409 (2014).CAS 
PubMed 
Article 

Google Scholar 
Simkanin, C. et al. Exploring potential establishment of marine rafting species after transoceanic long-distance dispersal. Glob. Ecol. Biogeogr. 28, 588–600 (2019).Article 

Google Scholar 
Raxworthy, C. J., Forstner, M. R. J. & Nussbaum, R. A. Chameleon radiation by oceanic dispersal. Nature 415, 784–787 (2002).CAS 
PubMed 
Article 

Google Scholar 
Stehli, F. G. & Webb, S. D. The great American biotic interchange., Vol. 4 (Springer Science & Business Media, 2013).Ronquist, F. Dispersal-vicariance analysis: A new approach to the quantification of historical biogeography. Syst. Biol. 46, 195–203 (1997).Article 

Google Scholar 
Ricklefs, R. E. & Bermingham, E. The concept of the taxon cycle in biogeography. Glob. Ecol. Biogeogr. 11, 353–361 (2002).Article 

Google Scholar 
Ma, H. An analysis of the equilibrium of migration models for biogeography-based optimization. Inf. Sci. 180, 3444–3464 (2010).Article 

Google Scholar 
Yiming, L., Niemelä, J. & Dianmo, L. Nested distribution of amphibians in the Zhoushan archipelago, China: can selective extinction cause nested subsets of species? Oecologia 113, 557–564 (1998).CAS 
PubMed 
Article 

Google Scholar 
Crisci, J. V., Katinas, L. & Posadas, P. Historical Biogeography: An Introduction. (Harvard University Press, 2003).Chen, R. et al. Adaptive innovation of green plants by horizontal gene transfer. Biotechnol. Adv. 46, 107671 (2021).CAS 
PubMed 
Article 

Google Scholar 
Schönknecht, G., Weber, A. P. & Lercher, M. J. Horizontal gene acquisitions by eukaryotes as drivers of adaptive evolution. BioEssays 36, 9–20 (2014).PubMed 
Article 
CAS 

Google Scholar 
Smith, A. B. Echinoderm phylogeny: morphology and molecules approach accord. Trends Ecol. Evolution 7, 224–229 (1992).CAS 
Article 

Google Scholar 
Bateman, R. M., Hilton, J. & Rudall, P. J. Morphological and molecular phylogenetic context of the angiosperms: contrasting the ‘top-down’ and ‘bottom-up’ approaches used to infer the likely characteristics of the first flowers. J. Exp. Bot. 57, 3471–3503 (2006).CAS 
PubMed 
Article 

Google Scholar 
Morris, J. L. et al. The timescale of early land plant evolution. Proc. Natl Acad. Sci. 115, E2274–E2283 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
Richter, S. The Tetraconata concept: hexapod-crustacean relationships and the phylogeny of Crustacea. Org. Diversity Evolution 2, 217–237 (2002).Article 

Google Scholar 
Dunn, C. W. et al. Broad phylogenomic sampling improves resolution of the animal tree of life. Nature 452, 745–749 (2008).CAS 
PubMed 
Article 

Google Scholar 
Caravas, J. & Friedrich, M. Of mites and millipedes: recent progress in resolving the base of the arthropod tree. BioEssays 32, 488–495 (2010).CAS 
PubMed 
Article 

Google Scholar 
Howard, R. J. et al. The Ediacaran origin of Ecdysozoa: integrating fossil and phylogenomic data. J. Geol. Soc. https://doi.org/10.1144/jgs2021-107 (2022).Newman, M. E. J. A model of mass extinction. J. Theor. Biol. 189, 235–252 (1997).CAS 
PubMed 
Article 

Google Scholar 
Cobbett, A., Wilkinson, M. & Wills, M. A. Fossils impact as hard as living taxa in parsimony analyses of morphology. Syst. Biol. 56, 753–766 (2007).PubMed 
Article 

Google Scholar 
Ruta, M., Krieger, J., Angielczyk, K. & Wills, M. A. The evolution of the tetrapod humerus: morphometrics, disparity, and evolutionary rates. Earth Environ. Sci. Trans. R. Soc. Edinb. 109, 351–369 (2018).
Google Scholar 
Puttick, M. N., Thomas, G. H. & Benton, M. J. High rates of evolution preceded the origins of birds. Evolution 68, 1497–1510 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
Sansom, R. S. & Wills, M. A. Fossilization causes organisms to appear erroneously primitive by distorting evolutionary trees. Sci. Rep. 3, 1–5 (2013).Article 

Google Scholar 
Brinkworth, A., Sansom, R. & Wills, M. A. Phylogenetic incongruence and homoplasy in the appendages and bodies of arthropods: why broad character sampling is best. Zool. J. Linn. Soc. 187, 100–116 (2019).Article 

Google Scholar 
Brown, J. W. & Smith, S. A. The past sure is tense: on interpreting phylogenetic divergence time estimates. Syst. Biol. 67, 340–353 (2018).PubMed 
Article 

Google Scholar 
Barba-Montoya, J., Dos Reis, M. & Yang, Z. H. Comparison of different strategies for using fossil calibrations to generate the time prior in Bayesian molecular clock dating. Mol. Phylogenetics Evolution 114, 386–400 (2017).CAS 
Article 

Google Scholar 
Sanderson, M. J. & Donoghue, M. J. Patterns of variation in levels of homoplasy. Evolution 43, 1781–1795 (1989).PubMed 
Article 

Google Scholar 
Alroy, J. Fossilworks: Gateway to the Paleobiology Database, http://fossilworks.org (2019).Benton, M. J. The Fossil Record 2. (Chapman & Hall, 1993).Cohen, K. M., Harper, D. A. T. & Gibbard, P. L. ICS International Chronostratigraphic Chart 2021/02, http://www.stratigraphy.org/ (2021).Gradstein, F. & Ogg, J. Geologic time scale 2004–why, how, and where next! Lethaia 37, 175–181 (2004).Article 

Google Scholar 
Rohde, R. A. The GeoWhen Database, (2005).O’Leary, M. A. et al. The placental mammal ancestor and the post–K-Pg radiation of placentals. Science 339, 662–667 (2013).PubMed 
Article 
CAS 

Google Scholar 
Kluge, A. G. A concern for evidence and a phylogenetic hypothesis of relationships among Epicrates (Boidae, Serpentes). Syst. Biol. 38, 7–25 (1989).Article 

Google Scholar 
Tolson, P. J. Phylogenetics of the boid snake genus Epicrates and Caribbean vicariance theory. Occasional Pap. Mus. Zool., Univ. Mich. 715, 1–68 (1987).
Google Scholar 
Clopper, C. J. & Pearson, E. S. The use of confidence or fiducial limits illustrated in the case of the binomial. Biometrika 26, 404–413 (1934).Article 

Google Scholar  More

Costanza, R. et al. Nature 387, 253–260 (1997).Article 

Google Scholar 
Toman, M. Ecol. Econom. 25, 57–60 (1998).Article 

Google Scholar 
Barbier, E. B., Burgess, J. C. & Folke, C. (eds) Paradise Lost? The Ecological Economics of Biodiversity (Routledge, 1994).
Google Scholar 
Daily, G. C. (ed.) Nature’s Services: Societal Dependence on Natural Ecosystems (Island, 1997).
Google Scholar 
Watson, R. T. et al. (eds) Global Biodiversity Assessment: Summary for Policy-Makers (Cambridge Univ. Press, 1995).
Google Scholar 
Reid, W. V. et al. Ecosystems and Human Well-Being: Synthesis (A Report of the Millennium Ecosystem Assessment) (Island, 2005).
Google Scholar 
Chichilnisky, G. & Heal, G. Nature 391, 629–630 (1998).Article 

Google Scholar 
Umaña Quesada, A. in Green Growth That Works: Natural Capital Policy and Finance Mechanisms from Around the World (eds Mandle, L. et al.) Ch. 13, 195–212 (Island, 2019).
Google Scholar 
Ouyang, Z. et al. in Green Growth That Works: Natural Capital Policy and Finance Mechanisms from Around the World (eds Mandle, L. et al.) Ch. 12, 177–194 (Island, 2019).
Google Scholar 
Salzman, J., Bennett, G., Carroll, N., Goldstein, A. & Jenkins, M. Nature Sustain. 1, 136–144 (2018).Article 

Google Scholar 
Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services. The Global Assessment Report on Biodiversity and Ecosystem Services: Summary for Policymakers (eds Díaz, S. et al.) (IPBES, 2019).
Google Scholar 
Chaplin-Kramer, R. et al. Science 366, 255–258 (2019).PubMed 
Article 

Google Scholar 
Hamel, P. et al. npj Nature Urban Sustain. 1, 25 (2021).Article 

Google Scholar 
Mandle, L., Ouyang, Z., Salzman, J. & Daily, G. C. (eds) Green Growth That Works: Natural Capital Policy and Finance Mechanisms from Around the World (Island, 2019).
Google Scholar 
Dasgupta, P. The Economics of Biodiversity: The Dasgupta Review. Abridged Version (HM Treasury, 2021).
Google Scholar 
Arkema, K. K. et al. Proc. Natl Acad. Sci. USA 112, 7390–7395 (2015).PubMed 
Article 

Google Scholar 
Ouyang, Z. et al. Proc. Natl Acad. Sci. USA 117, 14593–14601 (2020).PubMed 
Article 

Google Scholar  More

Illustration of the compound event indicesBuilding on earlier studies24,25, we develop two univariate indices to model concurrent climate conditions, i.e., a CWD index that varies from compound cold-wet conditions to CWD conditions, and a CCD index that varies from compound warm-wet conditions to CCD conditions (see “Methods”). The two indices incorporate the dependence between temperature and precipitation and are a measure of how warm/cold and dry a point is relative to the distribution of climate conditions at a given location. We illustrate the two indices on two grid points that have strong but opposite temperature-precipitation correlation. In the case where temperature and precipitation are strongly negatively correlated, the CWD index is well aligned with the primary axis of the bivariate distribution (Fig. 1a). In the case where temperature and precipitation are strongly positively correlated, the same holds for the CCD index (Fig. 1d). As illustrated for several concurrent hot-dry and cold-dry events that occurred around the globe, the two indices well capture these events (Supplementary Figs. 1 and 2).Fig. 1: The relationship between precipitation and temperature and compound indices.a Scatter plot of summer precipitation and temperature anomalies (z-score) with corresponding CWD index in color (see “Methods”). The location is at 97.25°W and 33.75°N from 1901 to 2018. b The same as a but for spring at 84.75°E and 66.75°N. c Same distribution as in a but colored based on the CCD index. d Same distribution as in b but colored based on the CCD index.Full size imageNotably, in the case where precipitation and temperature are strongly positively correlated, the CWD index indicates the relative anomalies of bivariate joint distribution, and some counterintuitive situations might occur relative to the univariate marginals (Fig. 1b). For instance, points might be labeled as strong CWD events (CWD index > 1.5) even though temperature is anomalously cold (temperature anomalies < 0, red dots in lower left quadrant of Fig. 1b). The CCD index exhibits similar behavior (Fig. 1c). This indicates an interesting property of the compound indices to identify strong compound conditions relative to bivariate distribution that are not necessarily extreme from a univariate perspective3,24,26,27.Widespread direct and lagged impacts of springtime compound climate conditionsTo evaluate the lagged summer vegetation responses to spring compound climate conditions, we compute partial correlation between CWD (CCD) spring and subsequent summer vegetation variation by controlling for the influence of summer compound climate conditions on these correlations (see “Methods”). Results show widespread negative associations between CWD spring and subsequent summer vegetation in the mid-latitudes (50°N).a–c The average standardized anomalies (z-score) of GPP during CWD spring but subsequent non-CWD summer (a), non-CWD spring but subsequent CWD summer (b), and consecutive CWD spring and summer (c) for areas in Fig. 2a where summer vegetation responds positively (r ≥ 0.22) to spring CWD climate conditions. d–f The same as a–c, but for soil moisture. g–i The same as a–c, but for runoff. The bar plots with dash lines (without dash line) indicate the average anomalies of multiple observation-based (model-based) products, and the circles indicate the average anomalies of each product. GLASS, LUE, NIRv, Flux-CRU, and Flux-ERA5 are observation-based GPP products, while model simulations are taken from TRENDYv6. GLEAM is observation-based soil moisture. GRUN represents observation-based runoff. GLDAS-VIC, GLDAS-Noah, GLDAS-Catchment, and FLDAS indicate assimilatory soil moisture and runoff that incorporate satellite- and ground-based observational products.Full size imageFig. 4: The responses of vegetation productivity and hydrological variables to CWD events in mid-latitudes (23.5–50°N/S).a–c The average standardized anomalies (z-score) of GPP during CWD spring but subsequent non-CWD summer (a), non-CWD spring but subsequent CWD summer (b), and consecutive CWD spring and summer (c) for areas in Fig. 2a where summer vegetation responds negatively (r ≤ −0.22) to spring CWD climate conditions. d–f The same as a–c, but for soil moisture. g–i The same as a–c, but for runoff. The bar plots with dash lines (without dash line) indicate the average anomalies of multiple observation-based (model-based) products, and the circles indicate the average anomalies of each product. For details on data see Fig. 3.Full size imageFig. 5: The effects of CCD events on vegetation productivity and hydrological variables in mid-to-high latitudes.a–c The average standardized anomalies (z-score) of GPP during CCD spring but subsequent non-CCD summer (a), non-CCD spring but subsequent CCD summer (b), and consecutive CCD spring and summer (c) for areas in Fig. 2b where summer vegetation responds negatively (r ≤ −0.22) to spring CCD climate conditions. d–f The same as a–c, but for soil moisture. g–i The same as a–c, but for runoff. The bar plots with dash lines (without dash line) indicate the average anomalies of multiple observation-based (model-based) products, and the circles indicate the average anomalies of each product. For details on data see Fig. 3.Full size imageCWD events increase vegetation productivity in high latitudesWe first analyze the direct responses of productivity to springtime and summertime CWD events across high latitudes ( >50°N, Fig. 3). Productivity increases during CWD spring and summer (Fig. 3a–c), which is consistent with vegetation responses (Supplementary Fig. 8a–c). Despite elevated spring greenness, spring water overall shows positive anomalies during CWD spring (Fig. 3d, f, g, i). This result indicates that spring greenness during CWD conditions is not associated with dry spring across high latitudes, which is further confirmed by similar anomalies in springtime TWS (Supplementary Fig. 8d, f). In contrast, severe water reduction is found in CWD summer (Fig. 3e, f, h, i). This suggests that despite the beneficial effects of CWD events on productivity in summer, they are associated with summer water deficit.Next, to analyze the lagged effects of springtime CWD events, we investigate the productivity anomalies in summer under three cases, namely CWD spring but non-CWD summer, non-CWD spring but CWD summer, and consecutive CWD spring and summer. Our results indicate that springtime CWD events have positive lagged effects on summer productivity across high latitudes (Fig. 3). Specifically, we find that during non-CWD summer (that is not favorable for summer vegetation growth) preceded by CWD spring, positive anomalies are still found in summer productivity (Fig. 3a). In contrast, during CWD summer (preceded by non-CWD spring), some models and observation-based products exhibit a reduction in summer productivity (Fig. 3b). We further find that summer productivity highly increases during consecutive events (Fig. 3c). Vegetation anomalies show similar behaviors (Supplementary Fig. 8a–c). Regarding the lagged responses of hydrological variables, CWD springs followed by non-CWD summers do not lead to water dryness, despite increased vegetation greenness (Fig. 3d, g). The magnitude of summer water deficit is similar for both cases that include CWD summer (Fig. 3e, f, h, i) and is consistent with summer TWS anomalies (Supplementary Fig. 8e, f). These results imply that in high latitudes, summer water reductions characterized by TWS, soil moisture, and runoff are not associated with increased spring greenness but are primarily caused by summer precipitation deficit.The productivity responses to compound climate conditions may be stronger than that to individual events through the synergistic effects of temperature and precipitation28. To investigate this, we compute the average anomalies in GPP and soil moisture associated with univariate events across the focus areas, which are then compared with the effects of CWD and CCD events in high latitudes (see “Methods”). Warm events can not only directly increase productivity but also show positive lagged effects (Supplementary Fig. 9a, b). In contrast, dry events reduce productivity (Supplementary Fig. 9e, f). This indicates that the direct and lagged positive effects of CWD events across high latitudes are mainly dominated by the warm component, while dry conditions have negative effects. Therefore, the warm-induced increase in productivity slightly exceeds that associated with CWD events (Supplementary Fig. 9b). Soil moisture under warm springs shows positive anomalies (Supplementary Fig. 9c, d), while they slightly decline during dry spring (Supplementary Fig. 9g, h). This suggests that the positive anomalies in soil water during CWD spring are driven by the warm component.CWD events reduce vegetation productivity in mid-latitudesHere, we first investigate the direct effects of springtime and summertime CWD events across mid-latitudes (23.5–50°N/S). Springtime productivity exhibits little changes during CWD spring (Fig. 4a, c), despite dry spring (Fig. 4d, f, g, i). When considering the direct effects of CWD events in summer, the results are similar, whereas the negative magnitude of productivity in summer is larger than that in spring (Fig. 4b, c). This difference suggests CWD conditions in summer show more adverse effects on productivity than that in spring in mid-latitudes. The anomalies in vegetation and TWS are consistent (Supplementary Fig. 10).Next, the lagged effects of springtime CWD events in mid-latitudes are assessed. In cases with CWD spring but non-CWD summer, summer productivity exhibits slight anomalies (Fig. 4a), with slightly decreased summer water (Fig. 4d, g). Summer productivity and water show much higher reductions in case with consecutive events (Fig. 4c, f, i) than for the case with only CWD summer (Fig. 4b, e, h). These results are supported by the responses of vegetation indices and TWS (Supplementary Fig. 10), revealing that springtime CWD events in mid-latitudes have negative lagged effects on summer productivity and water availability.The direct and lagged effects of individual events are finally compared with that of CWD events in mid-latitudes. Dry conditions in spring and summer directly decrease productivity and cause soil water dryness (Supplementary Fig. 11a–d). Moreover, dry spring depletes soil moisture earlier, which, in turn, causes dry summer and reduction in productivity during non-dry summer (Supplementary Fig. 11a, c). This indicates that dry springs have negative lagged effects on summer productivity. In contrast, productivity and soil water show positive anomalies during warm springs, while they show negative anomalies in summer (Supplementary Fig. 11e–h). These results suggest that the direct and lagged negative effects of CWD springs are dominated by the dry component in mid-latitudes, while the warm component mitigates the negative effects of the dry component in spring. Accordingly, the decline in productivity due to dry conditions thus exceeds that triggered by CWD events (Supplementary Fig. 11b).Decreased vegetation productivity due to the negative synergistic effects of CCD eventsHere, we first investigate the direct effects of CCD events across mid-to-high latitudes. Productivity reductions are found during springtime and summertime CCD events (Fig. 5a–c) concurrent with water reductions (Fig. 5). Vegetation and TWS show similar behaviors during CCD spring and summer (Supplementary Fig. 12). These results reveal that CCD events in spring and summer can impose direct adverse impacts on productivity and soil water across mid-to-high latitudes. The productivity reductions in spring and summer are similar in magnitude (Fig. 5a, b), indicating that CCD events between spring and summer can cause similar damage to productivity.We then analyze the lagged effects of springtime CCD events. Our results indicate that springtime CCD events show negative lagged effects on summer productivity and cause summer water reductions in mid-to-high latitudes (Fig. 5). Specifically, we find that in cases with CCD spring but non-CCD summer, summer productivity and water exhibit strongly negative anomalies (Fig. 5a, d, g). Moreover, summer anomalies are higher during consecutive events (Fig. 5c, f, i) than the cases including only CCD summer (Fig. 5b, e, h). Vegetation indices and TWS show similar responses (Supplementary Fig. 12). Our results further indicate that CCD spring has more severe negative lagged effects on productivity than CWD spring. That is, we find that in comparison to cases with preceding CWD spring and consecutive CWD events, summer productivity shows higher reduction in cases with preceding CCD spring and consecutive CCD events (Fig. 4a, c versus Fig. 5a, c). Moreover, in cases with CCD spring but non-CCD summer (Fig. 5a, d, g), summer anomalies are close to those in scenarios with non-CCD spring but CCD summer (Fig. 5b, e, h). The vegetation and TWS anomalies further confirm this situation (Supplementary Fig. 12a, b, d, e). These results suggest that the lagged effects of CCD spring can be of similar magnitude as their direct adverse effects.We finally compare the direct and lagged effects of individual events with that of CCD events in mid-to-high latitudes. Cold conditions in spring and summer directly reduce productivity but show weak effects on soil moisture (Supplementary Fig. 13a–d), and cold spring shows negative lagged effects on summer productivity (Supplementary Fig. 13a). Dry events show direct and lagged negative effects on productivity and soil moisture (Supplementary Fig. 13e–h). These results imply that the negative lagged effects of CCD springs are dominated by both cold and dry components. The effects of CCD events on productivity mostly exceeds the individual dry or cold impacts (Supplementary Fig. 13a, b, e, f). More

OverviewIn this study, we used VIIRS active fire detections to track the dynamic evolution of all fires in California from 2012 to 2020 (Fig. 1). We developed an approach that has the following steps. First, after reading the satellite fire pixel data at each 12-hour time step, the new fire pixels are aggregated into multiple clusters using the fire pixel locations and an automatic clustering algorithm. These clusters are then spatially compared to existing fire objects. If a cluster is not close to any existing active fire object, we use all fire pixels within the cluster to form a new fire object. If a cluster is located near an existing fire object which is still active, we view the cluster as an extension of the existing fire. In this case, we append all pixels within the cluster to the corresponding existing fire object, allowing the existing object to grow. When a fire expands and gets close enough (within a pre-defined distance threshold) to an existing active fire object, we merge the two objects. For each time step (12 hours in this case for the two overpasses), we derive or update a suite of attributes and status indicators associated with each fire event, including pixel-level attributes of fire and surface properties, vector geometries related to the fire shape, and meta-attributes characterizing the entire fire object.Data inputSatellite remote sensing instruments provide active fire detections with accurate geographical location and broad spatial coverage. The primary data for this fire tracking system are active fire locations and the fire radiative power (FRP) recorded by the VIIRS instrument aboard the Suomi National Polar-orbiting Partnership (Suomi-NPP) satellite24. VIIRS observes Earth’s surface twice each day in low and mid latitude regions, with local overpass times of approximately 1:30 am and 1:30 pm. Compared to its predecessor, the MODIS sensors on the Terra and Aqua satellites, VIIRS has a higher spatial resolution and can detect smaller and cooler fires24. Also, the VIIRS instrument provides a more consistent pixel area across the image swath25, resulting in more accurate estimates of active fire location. Therefore, compared with MODIS, the VIIRS active fire products can be used to map fire event progression with higher accuracy21. Two streams of VIIRS active fire data are operationally produced using a contextual fire detection algorithm24, drawing upon VIIRS moderate resolution band (M-band) and imaging band (I-band) reflectance and radiance data layers. In this fire tracking system, we used the Suomi-NPP VIIRS I-band fire location data product (VNP14IMGML, Collection 1 Version 4) that contains the centre location, FRP, scan angle, and other attribute fields associated with each pixel. The I-band fire detection product has a 375-m spatial resolution at nadir (the sub-satellite point) and an average resolution across the full swath of about 470 m. Theoretical estimates of fire detection efficiency for the VIIRS sensor indicate that during the day, VIIRS can detect 700 K fires with 50% probability that have a size of about 200 m2 (a 15 m × 15 m fire area)24. During night, the detection efficiency increases, and VIIRS can detect 700 K fires as small as 40 m2. From a fire spread tracking perspective, these detection efficiencies imply that in many instances, the area of a fire pixel that is covered with flaming fire combustion is several orders of magnitude smaller than the overall pixel size. The VNP14IMGML data, available from 2012 onwards, were downloaded from the University of Maryland VIIRS Active Fire website (https://viirsfire.geog.umd.edu/).Land cover data are an additional input in the system required to classify different fire types and determine the spatial connectivity threshold. Here we use the U.S. National Land Cover Database (NLCD 2016)26 that is available from the Multi-Resolution Land Characteristics (MRLC) Consortium website (https://www.mrlc.gov/national-land-cover-database-nlcd-2016). We aggregated the original 30-m data to match the spatial resolution of VIIRS active fire data, and merged the original 16 classes into several groups: ‘Water’, ‘Urban’, ‘Barren’, ‘Forest’, ‘Shrub’, ‘Grassland’, and ‘Agriculture’. We used the 1000-hour dead fuel moisture from the high-resolution (4 km) gridMET product27 for the purpose of separating wildfires and management fires. This gridMET dataset was computed from 7–day average conditions composed of day length, hours of rain, and daily temperature and humidity ranges. Regularly updated gridMET data are available from the Climatology Lab website (http://www.climatologylab.org/gridmet.html).Other ancillary and validation datasets used in this study included a shapefile of California borders and fire perimeters from the California Forestry and Fire Protection’s Fire and Resource Assessment Program (FRAP) dataset (https://frap.fire.ca.gov/mapping/maps/).Fire object hierarchyFire detections from VIIRS are dynamically tracked within the framework of a three-level object hierarchy (Fig. 1). The lowest level is the fire pixel object, which includes the geographical location (latitude and longitude), the FRP value, and the origin (first assigned fire object id). The second level is the fire object, which includes all attributes associated with each individual fire event at a particular time step (Table 2). Each fire object includes one or more fire pixel objects, a unique identification number (id), and a set of attributes associated with the whole fire. Two types of fire attributes are derived and recorded for each fire object. The first type encompasses temporal (e.g., ignition time, duration) and spatial (e.g., centroid, ignition location) characteristics of the object as well as general properties (e.g., size, type, active status). The second type is the geometric information related to the fire object, including the fire perimeter, the active fire front line, and the newly detected fire pixel locations (stored as vectors). All fire objects in the State of California are combined to form an allfires object, to characterize the whole-region fire situation at a specific time step. The allfires object comprises a list of fire objects, and also contains meta information representing the statistics of all fires and the records describing fire evolution. A full list of the attributes associated with the pixel object, the fire object, and the allfires object is presented in Table 2.Table 2 List of main attributes associated with pixel, fire and allfires objects.Full size tableFire event trackingThe fire records (locations and FRPs) from the monthly VIIRS active fire location products (VNP14IMGML) are read into the system at each half-daily time step (roughly 1:30 am and 1:30 pm local time). We apply spatial and temporal filters to the data to extract active fire pixels recorded in California during each 12-hour time interval. We also apply quality flag filters (thermal anomaly type of ‘0: presumed vegetation fire’ in VNP14IMGML)) to ensure the use of only pixels likely associated with vegetation fires. The fire location and FRP values are used to create fire pixel objects. To speed up the calculation, the newly detected active fire pixels after filtering are first aggregated to specific clusters using the distances between them and an automatic clustering algorithm. In this initial aggregation algorithm, a ball tree28 is created to partition all newly detected active fire pixels into a nested set of hyperspheres in a 2-D space (latitude and longitude). This space partitioning data structure can be used to expedite nearest neighbours search29 and allow for quick cluster grouping. Here we refer to a cluster as a collection of pixel objects that are recorded at the same time step and are also spatially nearby. In the following steps, all pixels within a cluster are considered as a whole for fire merging and creation.We define an extended area for every existing fire object as the fire vector perimeter (see the section of Calculating and recording fire attributes for detail) plus a radial buffer that depends on the fire type property of the object. The buffer is set to 5 km for forest fires and 1 km for other fire types (shrub, crop, urban), considering that the fire spread rate can differ across biomes13. We then evaluate the spatial distance between the perimeters of a newly classified cluster and all existing active fire objects (a fire object keeps an active status if one or more active fire pixels associated with it are detected during the past 5 days), and calculate the shortest distance. If the shortest distance is smaller than the buffer of the associated existing active fire (i.e., new cluster overlaps with the extended area of an existing fire object), we assume all fire pixels in the new cluster are associated with the growth of the existing fire object at the current time step (Fig. 2). The existing fire object is updated by appending all fire pixel objects within the new cluster. If a newly classified cluster does not overlap with the extended area of any existing active fire object, we assume this is a new fire. A new fire object (by assigning a new fire id) is created using all fire pixel objects in the cluster.With the addition of new fire pixels, an existing fire object may expand and touch the extended area of another existing active fire object. If this happens, we assume that these two existing fire objects merge into a single object at this time step. All fire pixels in the fire object with a higher id number (a later start date, termed as the ‘source fire’) are appended to the fire object with lower id number (earlier start date, termed as the ‘target fire’) in this case. We record the id of the target fire in a list of fire mergers, and update all attributes associated with this fire (Fig. 3). In order to avoid double counting, the source fire object (with all pixels being transferred to the target fire object) is flagged as invalid, and is excluded from statistical analysis of fire events.Fig. 3The time series of growth for the SCU Lightning Complex fire (2020). Panel (b) shows the fire size of the SCU fire (total area within the fire object perimeter) at half-daily time steps. A fraction of the fire growth (shown in orange) was due to the addition of newly detected fire pixels. Panel (a) shows the number of new fire pixels (associated with the SCU fire object) detected at each time step. The other part of the fire growth (shown in red) was due to the merging with existing fire objects. Panel (c) shows the number of fire pixels in the existing objects that were merged to the SCU fire object.Full size imageCalculating and recording fire attributesOther than individual fire pixels contained in a fire object, several core attributes (properties and geometries) are also dynamically updated at each time step and are used for fire tracking and characterization.Important time-related attributes include the fire ignition time (the time step at which the first fire pixel within the fire object was detected), the fire end time (the latest time step with an active fire observation), and the fire duration (the time difference between the ignition time and end time). If a fire object does not have new active fire pixels appended during 5 consecutive days (i.e., the fire end time is more than 5 days before the present time step), its status is set to inactive. Once inactive, a fire object is no longer evaluated for use in future clustering (i.e., new active fire detections later will form new fire objects, even if they are spatially close to the inactive fire object).Each fire object is assigned to a specific fire type. The fire type is identified using the major land cover type within the fire perimeter (Table 3). In an initial analysis, we found that prescribed fires, on average, have higher coarse fuel moisture levels than wildfires. Therefore, we also record the 1000-hour fuel moisture (fm1000) from the gridMET dataset27 for each fire object (corresponding to the ignition time step) and use this value to divide forest and shrub fires further to wildfire and prescribed types.Table 3 Classification of fire types based on dominant land cover type (from the US National Land Cover Database) within each fire perimeter and the 1000-hr fuel moisture (FM-1000, from gridMET dataset) at the time of ignition.Full size tableAn essential step in this object-based fire tracking system is to determine the vector shape of the fire perimeter. In this system, we use an alpha shape30 algorithm to derive bounding polygons containing fire pixels in a fire object. For an alpha shape, the radius of the disks forming the curves in the polygon is determined by the alpha parameter α. Compared with the commonly used convex hull, the alpha shape hull is able to capture the irregular shapes around the fire perimeter more accurately22.To identify the optimal values for the α parameter, we performed the following analysis. First, we derived the final fire perimeters for all large fires that occurred in California during the 2018 wildfire season using a set of α values ranging from 500 m to 10 km and compared the results with more refined fire perimeters from the Fire and Resource Assessment Program (FRAP) dataset (Fig. 4). Large magnitude α values tended to overestimate the total burned area, while small α values often fragmented a large fire event. We found that a value of α = 1 km was optimal in terms of balancing the ability of the hull to catch the boundary shape and to keep the integrity of a fire object. For each time step, we applied the alpha shape algorithm to all fire pixel locations associated with a fire object since the time of ignition. This processing step resulted in a concave hull with the shape of polygon or multipolygon. To account for the pixel size, we expanded the concave hull to the fire perimeter using a buffer size equal to half of the VIIRS nadir cross-track pixel width (187.5 m). The alpha shape algorithm does not work when the total number of fire pixels (npix) is less than 4. If npix equals 3, we used a convex hull algorithm and the same 187.5 m buffer to determine a polygon perimeter. If npix is 1 or 2, circles centered on the fire pixel location with radius of 187.5 m were used.Fig. 4Optimization of the alpha shape parameter (α). For all large fires (final size  > 4 km2) in California during 2018, fire perimeters were estimated using VIIRS active fires and different alpha parameters. By comparing (a) the burned area (BA) and (b) the number of fire objects with the FRAP data, an optimal alpha parameter of 1 km was identified for use in this study (shown in red). The vertical bars and lines show the mean and 1-std variability from all fires. The dashed blue lines indicate the ideal values when compared to FRAP. Panels (c)–(h) show the fire perimeters derived using different alpha shape parameters for two sample fires in 2018. The shapes with pink color are final FEDS fire perimeters derived from VIIRS active fires using the alpha shape algorithm. The blue shapes represent the corresponding fire perimeters from the FRAP dataset. Overlap between FRAP and FEDS is shown in purple.Full size imageWe also calculate the active front line for each fire object at each time step. The active fire front consists of the segments of the fire perimeter that are actively burning and releasing energy and emissions. The position of the active fire line is critical in evaluating the fire risk, estimating the fire emissions, and predicting fire spread. We derive the active portion of the fire perimeter as segments that are within a 500 m radius of newly detected fire pixel locations. We found that this threshold allowed for a continuous projection of the active fire front in rapidly expanding areas of large wildfires during the 2018 fire season; this threshold may be optimized in future work to maximize performance metrics for fire model forecasts. The resulting active line for each fire at each time step has the shape of a linestring (object representing a sequence of points and the line segments connecting them), a multi-linestring (a collection of multiple linestrings), or a linear ring (closed linestring). Figure 5 shows an example map of the fire perimeters and active fire front lines on September 8 during the 2020 wildfire season.Fig. 5An example map of fire perimeters and fire active fronts in California. The map was created using the fire event data suite (FEDS) as of the Suomi-NPP afternoon overpass (~1:30 pm local time) of Sep 8, 2020. The background is the Aqua MODIS Corrected Reflectance Imagery (true color) recorded at the same day (provided by the NASA Global Imagery Browse Services). The active front line of a fire is shown in yellow, active fire areas are shown in red, and the area of inactive (extinguished) fires are shown in dark red.Full size imageAdditional fire properties, such as the fire area and active fire line length, are also derived using these geometries of the fire object (see Table 2). Note this list can be easily expanded to include more user-defined properties with the help of the fire object core vector data.The allfires object contains a list of all existing fire objects at a time step. This object also records the ids of fire objects that have been modified (including fires newly formed, fires that expanded with new pixel additions, fires with pixels addition due to merging, and fires that just became invalid) at the current time step.Creating the fire event data suite (FEDS)By tracking the spatiotemporal evolution of all fire objects in California, we derived a complete dataset of fire events for each calendar year (Jan 1 am – Dec 31 pm) during the Suomi-NPP VIIRS era (2012–2020). The dataset contains four products that represent the fire information in California at multiple spatial scales and from different perspectives (Fig. 1 and Table 4), ranging from the most detailed and memory-intensive data format (Pickle) to the most high-level format (CSV).Table 4 Data structure of the FEDS.Full size tableThe first product is the direct serialization result of the allfires object at each time step (twice per day). The product is stored as a Pickle file31 which allows for analysis of the complex allfires object structure (including all attributes associated with all fire objects it contains). This file also serves as the restart file for continued fire tracking at any time step, which is essential for the operational mode using the near-real-time fire data. By restoring an exact copy of the previously pickled allfires object, any attribute in the allfires object can be deserialized from the saved files. The Pickle file is the most basic data product in the dataset, and is created at each half-day time step.The second product (Snapshot) represents a more accessible and self-explanatory variant of the Pickle serialization product. In this product, we tabulated important diagnostic attributes for each fire and saved them in GeoPackage32 data files. Each GeoPackage file includes three data layers: one contains the properties and the fire perimeter geometry, another contains the active fire line geometry, and a third contains the new fire pixel location geometry. This product, created at a half-daily time step, allows for a more straightforward interpretation of regional fire status at a particular time step. We also created a GeoPackage file that summarizes the final fire perimeters and attributes for all fires during the whole study period (2012–2020).The third product (Largefire) focuses on the temporal evolution of individual large fires with an area greater than 4 km2. At each time step, the time series of properties and geometries (fire perimeter, active fire line, and new fire pixel locations) for each of the large fires are extracted and saved to GeoPackage files. This product facilitates the visualization and analysis for an individual targeted fire (Fig. 6) and is particularly useful in the near-real-time evaluation, forecasting, and policy making.Fig. 6The spatiotemporal evolution of the Creek fire (2020). Contours and dots reflect the fire perimeters and newly detected fire pixels at each 12-hour time step. Data for the period of Sep 5 am–Nov 6 am, 2020 are shown.Full size imageThe fourth product (Summary), which is stored as NetCDF and CSV files and created at the end of a fire season, records the all-year time series of fire statistics (including major fire attributes such as number, size, duration, fire line length, etc.) over the whole State of California. This product provides a feasible regional summary of the temporal evolution of fires.Potential for near-real-time (NRT) fire event trackingWhile the main objective of this paper is to apply the object-based fire tracking system to historical VIIRS fire detections and create a retrospective multi-year FEDS, we note that this system has the potential to be used for tracking fire events in near-real-time, providing rich and valuable information for fire management and short-term risk assessment. We have experimented with the use of this system for NRT fire event tracking in California using the daily NRT Suomi-NPP VIIRS active fire detection product (VNP14IMGTDL, collection 6) as the main data source. The VNP14IMGTDL product is routinely produced and is publicly available at the NASA Fire Information for Resource Management System (FIRMS). Since the NRT product undergoes less rigorous quality assurance, we use only fires with ‘nominal’ or ‘high’ confidence levels from the NRT product for fire tracking. Some active fire detections from the NRT data are potentially associated with static non-vegetation fires (e.g., fires from gas flaring in oil and gas or landfill industries or false detections due to reflection from solar panels) and are not the main interest for vegetation fire studies. To avoid the unnecessary computation associated with these static fires, we record and evaluate the fire pixel density for each fire object at each time step. When a small fire ( 20 per km2), it is considered to be a static fire and subsequently labelled as invalid.Similar to the retrospective FEDS, we use the active fire detections to create an object serialization product, a regional snapshot GIS product, and a time series product of large fire evolution twice daily. This experimental NRT data will be available upon publication through a university hosted server. More

Study areaGrassland herbage samples are from Shaerqin base, institute of grassland research of CAAS (Chinese Academy of Agricultural Sciences). We obtained the permission of the institution to take HSI of the grassland sample. Our work did not cause damage to grassland. Researcher Weihong Yan of the institute provided us with relevant information about grassland. The land use type in the study area is mainly grassland, which is composed of forage species, most of which are representative species of typical grassland. We take this area as an example to conduct research on grass classification. By enriching the relevant recognition technology, it can also be used as a reference for the pastures of other grasslands. The grass species Grass1 for the experiment is shown in Table 1. The official introduction of plant materials is detailed in the flora of China15.Table 1 Samples information for Grass1 dataset.Full size tableThe field hyperspectral platformWe assemble a system for collecting HSI in the field: HyperSpec©PTU-D48E HSI instrument, high-precision scanning PTZ, tripod, data analysis software Hyperspec, etc. The light source is natural light. The imaging instrument is in line scanning mode. Table 2 shows the technical parameters.Table 2 Technical parameters of hyperspectral instrument.Full size tableData collectionIn July 2021, the data was collected during the lush grass growth period. Collect data from 11:00 a.m. to 2:00 p.m. every day. At this time, it is sunny, cloudless and the wind force does not exceed level 2. So as to ensure the consistency of the acquisition time line and avoid the influence of different degrees of light on the reflectivity as far as possible. The measuring points are arranged facing the sun and the opposite direction of the shadow. We collect data from different angles of the grassland, which is based on the growth of various types of forages, and selects relatively concentrated places within the study area. Each shot is a single category of grass. The image resolution is 1166 × 1004 pixels (Fig. 1). The imaging spectrometer is fixed with scanning head when shooting. Data acquisition and transmission are executed on Hyperspec software. Then save it as a BIL file. The ENVI5.3 software was used to extract the forage spectrum to establish the dataset Grass1. Well balanced regions with a clear image, uniform spectral distribution are selected for further segmentation. The average value of spectral reflectance of grass pixels was taken as the reflectance spectrum of a single type of grass.Figure 1True color map of grass samples.Full size imageMethodologyIn Fig. 2, we present the framework of visible-NIR hyperspectral classification of grass based on multivariate smooth mapping and extreme active learning (MSM–EAL). Specifically, we first introduce the proposed MSM algorithm for global enhanced spectral reconstruction, which utilizes smooth manifold projection technology to alleviate the problems of difficult feature selection and redundant data. Then, the EAL framework is proposed to address the matter of hyperspectral labeled samples and spectral classification. In the following, each step of this method will be presented in detail.Figure 2Proposed MSM–EAL framework for grass HSI classification.Full size imageThe proposed MSM algorithmIn the process of field HSI acquisition, on the one hand, the surface distribution of grass is uneven and the plant height is different, causing certain scattering effect and coverage spectrum change. On the other hand, HSI is easy to be disturbed by external natural factors such as light, wind and shadow, resulting in a certain degree of distortion. Multiplicative scatter correction (MSC) is a scattering correction effect, which helps to eliminate the scattering effect caused by the above reasons and enhance the spectral variability. The moving window smooth spectral matrix (Nirmaf) belongs to the smooth effect, which improve the signal-to-noise ratio of the spectrum and reduce the influence of random noise16,17. Preprocessing methods are different and related to each other. We design an enhanced preprocessing multivariate smooth (MS) method that fusing MSC and smooth Nirmaf to target grass spectral signal features. In the follow-up, a model will be established to verify the validity of MS.Most of the high-dimensional spatial data have the characteristics of being embedded in a manifold body, so the manifold learning isometric feature mapping (Isomap) based on spectral theory is adopted. Isomap preserves the global geometric features of the initial data and extracts features by reconstructing the underlying smooth manifold of HSI. It is nonlinear dimensionality reduction based on linear and multidimensional scaling transformation18. Isomap has been applied in image and HSI classification19,20, but there is no report on visible-NIR hyperspectral classification of grass.In view of the above, we proposed the multivariate smooth mapping (MSM) spectral reconstruction algorithm, which can be represented as follows:$$ MSM_{z} { } = { }frac{{left( {P_{j} – b_{j} } right)left( {2n + 1} right) + n_{j} cdot mathop sum nolimits_{j = – n}^{n} C_{j} P_{k + j} }}{{n_{j} left( {2n + 1} right)}} + V_{Z} F_{Z}^{frac{1}{2}} { } $$
(1)
where Pj, bj, and Cj represent the raw reflectance value of spectrum j, baseline shift amount, and weight factor, respectively, k and nj represent the polynomial degree and offset, respectively. MSMz is the feature cube reconstructed to Z dimension from the spectrum calculated by 2n + 1 moving window width, V eigenvector matrix and F eigenvalue matrix.In Isomap equidistant mapping, the shortest path of edge Pi Pj needs to be solved, and the representation matrix is:$$ D_{G} = [d_{G}^{2} (P_{i} ,P_{j} )]_{i,j = 1}^{n} $$
(2)
where d (Pi, Pj) is the weight of the edge Pi Pj calculated from the neighborhood graph G and its side Pi Pj.The proposed EAL frameworkLabeling hyperspectral samples is expensive in terms of time and cost, at the same time, the lower spatial resolution and more bands increase the difficulty of labeling. Active learning (AL) provides an efficient labeling strategy, which only needs to label a relatively small number of samples to learn a more accurate model21. The pool-based AL selects the most informative samples according to the query strategy for limited labeling through iteration, so as to facilitate model improvement. Commonly used query strategies are uncertainty criteria, such as least confidence22, the bayesian active learning disagreement (BALD), the entropy sampling23, etc.Due to there is still an over-fitting problem, different strategies such as hybrid prediction and regularization need to be used for non-recursive datasets24. The research25 proposed that extreme gradient boosting algorithm (XGBoost) based on gradient boosting. As a classification method, XGBoost has been successfully applied in Kaggle competition and other fields. Its most important feature for visible-NIR hyperspectral classification is that can easily and directly classify according to features, and the physical interpretation of features can help understand the electronic nature behind spectral classification. XGBoost is a machine learning algorithm based tree structure that integrates multiple weak classifiers to achieve flexible and high-precision classification. It is an upgraded version of gradient boosting decision tree. The optimization process of XGBoost entailed: (1) Expanding the objective function to the second order, and finds a new objective function for the new base model to improve the calculation accuracy. (2) L2 regularization term is added to the loss function to prevent over-fitting. (3) Using blocks storage structure realize automatic parallel computing26,27. The algorithm steps are as follows:The objective function:$$ Lleft( Phi right) = mathop sum limits_{i} lleft( {y^{i} ,widehat{{y^{i} }}} right) + mathop sum limits_{k} Omega left( {f_{k} } right) $$
(3)
In formula (3), the first and second terms are the loss function term and the regularization term, respectively. Where,$$ Omega left( {f_{k} } right) =upgamma {text{T}} + frac{1}{2}lambda left| w right|^{2} $$
(4)
γ and λ are regularization parameters which are used to adjust complexity of the tree.Next, second derivative Taylor expansion of the objective function. Where (g_{i}) and (h_{i}) are the first derivative and second derivative, respectively.$$ L^{left( t right)} = mathop sum limits_{i = 1}^{n} lleft( {y_{i} ,widehat{{y_{i}^{t – 1} }} + f_{t} left( {x_{i} } right)} right) + Omega left( {f_{t} } right) $$
(5)
$$ g_{i} = partial_{{hat{y}_{i} (t – 1)}} lleft( {y_{i} ,widehat{{y_{i}^{t – 1} }}} right) $$
(6)
$$ h_{i} = partial_{{widehat{{y_{i} }}(t – 1)}}^{2} lleft( {y_{i} ,widehat{{y_{i}^{t – 1} }}} right) $$
(7)
$$ {text{L}}^{left( t right)} approx mathop sum limits_{i = 1}^{n} left[ {lleft( {y_{i} ,widehat{{y_{i}^{t – 1} }}} right) + g_{i} f_{i} left( {x_{i} } right) + frac{1}{2}h_{i} f_{t}^{2} left( {x_{i} } right)} right] + Omega left( {f_{t} } right) $$
(8)
Final objective function:$$ {hat{text{L}}}^{ i} left( q right) = – frac{1}{2}mathop sum limits_{j = 1}^{T} frac{{(mathop sum nolimits_{{i in I_{j} }} g_{i} )^{2} }}{{mathop sum nolimits_{{i in I_{j} }} h_{i} + lambda }} + gamma T $$
(9)
Equation (9) can be used as the fraction of tree cotyledons, and the tree structure is directly proportional to the fraction. If the result after splitting is less than the maximum value of the given parameter, the cotyledon depth stops growing24,28.AL solves the problems of limited number and high cost of grass hyperspectral labeling samples. The default model of traditional AL is logistic regression, which is mostly studied on the ideal public dataset. However, the actual data has more uncertain noise, which still poses a certain challenge to AL. Consequently, we propose the extreme active learning (EAL) framework to minimize the classification cost of visible-NIR hyperspectral. The framework replaces the logistic regression model with XGBoost. Taking advantage of AL, XGBoost can improve performance with less training marker samples. By jointing of XGBoost and AL, EAL provides significantly better results than AL in field Grassl dataset recognition. Additionally, based on the characteristics of XGBoost, EAL more intuitively enhances the physical essence behind spectral classification than AL. Algorithm 1 summarizes the workflow of EAL framework.Random forest (RF) and decision tree (DT) were used to compare with EAL. RF and DT are frequently used in the field of grassland remote sensing9,29. Furthermore, RF, DT and XGBoost have the same point is that are learning algorithms based on tree structure. DT determines the direction by judging the conditions of the decision node12. RF is an integrated learning of multiple decision trees30. More