# Top 10+1 indicators for assessing forest ecosystem conditions: A five-decade fragmentation analysis

Bruna Almeida, Pedro Cabral, Catarina Fonseca, Artur Gil, Pierre Scemama

<figure><img src="image_1.png" /></figure>

## To cite this version:

Bruna Almeida, Pedro Cabral, Catarina Fonseca, Artur Gil, Pierre Scemama. Top 10+1 indicators for assessing forest ecosystem conditions: A five-decade fragmentation analysis. Science of the Total Environment, 2024, 957, pp.177527. 10.1016/j.scitotenv.2024.177527. hal-05244728

HAL Id: hal-05244728

https://hal.univ-brest.fr/hal-05244728v1

Submitted on 8 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

# Top 10+1 indicators for assessing forest ecosystem conditions: A five-decade fragmentation analysis

Bruna Almeida<sup>b, c, *</sup>, Pedro Cabral<sup>a, b, *</sup>, Catarina Fonseca<sup>c, d</sup>, Artur Gil<sup>c, e</sup>, Pierre Scemama<sup>f</sup>

## HIGHLIGHTS

* Cellular automata models large-scale forest resources and predict future scenarios

* The first principal component captures the primary trend describing forest dynamics

* Indicators reflect landscape structure, fragmentation, connectivity and complexity

* Fragmentation leads to more patches, smaller patch sizes, and reduced connectivity

* Findings have implications for landscape development and resource management

## GRAPHICAL ABSTRACT

<figure><img src="image_5.png" /></figure>

## ARTICLE INFO

Editor: Manuel Esteban Lucas-Borja

**Keywords:**

Forestry
Ecological applications
Environmental conservation
Corine land use land cover
Carbon storage and sequestration
GIS modelling

## ABSTRACT

Globally, land use change has consistently resulted in greater losses than gains in aboveground biomass (AGB). Forest fragmentation is a primary driver of biodiversity loss and the depletion of natural capital. Measuring landscape characteristics and analyzing changes in forest landscape patterns are essential for accounting for the contributions of forest ecosystems to the economy and human well-being. This study predicts national forest distribution for 2036 and 2054 using a Cellular Automata (CA) system and assesses ecosystem conditions through landscape metrics at the patch, class, and landscape levels. We calculated 130 metrics and applied a Variance Threshold method to remove features with low variance, testing different thresholds. The first filtered-out metrics were further analysed through Principal Component Analysis combined with a Feature Importance technique to select and rank the top 10 indicators: effective mesh size, splitting index, mean radius of gyration, largest patch index, mean core area, core area percentage, Simpson's evenness index, mutual information, Simpson's diversity index, and mean contiguity index. The eleventh selected indicator is the AGB density, a structural measurement for ecosystem condition and a proxy for forest carbon storage and sequestration assessments. From 2000 to 2018, the national AGB forest carbon stock decreased from 131.5 to 91.3 Megatons (Mt)

<footnote><sup>a</sup> School of Remote Sensing and Geomatics Engineering, Nanjing University of Information Science and Technology, Nanjing 210044, China</footnote>

<footnote><sup>b</sup> NOVA Information Management School (NOVA IMS), Universidade Nova de Lisboa, Campus de Campolide, 1070-312 Lisboa, Portugal</footnote>

<footnote><sup>c</sup> Centre for Ecology, Evolution and Environmental Changes (cE3c); Azorean Biodiversity Group (GBA), University of the Azores (UAc), Faculty of Sciences and Technology (FCT-UAc); Rua Mãe de Deus, Campus Universitário de Ponta Delgada, Edifício do Complexo Científico, 9500-321 Ponta Delgada, Portugal</footnote>

<footnote><sup>d</sup> MARE - Marine and Environmental Sciences Centre/ARNET - Aquatic Research Network, Faculdade de Ciências, Universidade de Lisboa, Portugal</footnote>

<footnote><sup>e</sup> IVAR - Research Institute for Volcanology and Risk Assessment, University of the Azores, 9500-321 Ponta Delgada, Portugal</footnote>

<footnote><sup>f</sup> IFREMER, Univ Brest, CNRS, UMR AMURE, Rue Dumont d'Urville, 29280 Plouzané, France</footnote>

<footnote>* Corresponding authors at: NOVA Information Management School (NOVA IMS), Universidade Nova de Lisboa, Campus de Campolide, 1070-312 Lisboa, Portugal. E-mail addresses: balmeida@novaims.unl.pt, cabral@nuist.edu.cn (B. Almeida).</footnote>
with expected values for 2036 and 2054 being 71.8 and 55.3 Mt., respectively. Landscape measurements quantitatively describe forest dynamics, providing insights into the structure, configuration, and changes characterizing landscape evolution. This research underscores the capability of CA models to map large-scale forest resources and predict future development scenarios, offering useful information for conservation and environmental management decisions. Additionally, it provides measurements to support Ecosystem Accounting by assessing forest extent and indicators of its conditions.

# 1. Introduction

Growing concerns about the degradation of natural capital, coupled with the increasing demand for ecosystem services (ES) underscore its importance in policy and governance frameworks (Farrell et al., 2021). Natural capital accounting stands out as a potential key tool to facilitate a shared understanding among stakeholders regarding resource management, environmental impacts, and the sustainability of land use decisions (Fleming et al., 2022). The adoption of the System of Environmental-Economic Accounting-Ecosystem Accounting (SEEA-EA) by the United Nations Statistical Commission in March 2021, marked a significant milestone in explicitly acknowledging how natural resources are accessed, utilized, and protected (Edens et al., 2022). This system moves beyond traditional Gross Domestic Product (GDP) accounting by integrating biodiversity and ecosystems into national economic planning and governance structures (Bateman and Mace, 2020). The SEEA-EA has the potential to strengthen laws, regulations, and policies, thereby promoting sustainable management of natural resources and balancing ecological integrity with human needs (King et al., 2024). To achieve this, the accounts involve (i) identifying ecosystem assets, (ii) measuring physical ecosystem accounts, such as the extent and indicators of its condition, and (iii) constructing monetary accounts including the valuation of ES (Bordt, 2018; Hein et al., 2020). Assessing ecosystems' capacity to provide goods and services in biophysical units is therefore central and is typically conducted using spatial models (Rijal et al., 2021).

Despite current ongoing efforts fostered by the United Nations Kunming-Montreal Global Biodiversity Framework, whose targets were recently adopted at the COP 15 to the Convention on Biological Diversity (CBD-UN) and are directly linked to the ES accounting concepts of the SEEA-EA, natural capital stocks and biodiversity have been deteriorated worldwide at unprecedented rates, demanding more effective actions (King et al., 2024). Thus, one of the goals agreed among the parties was to maintain, enhance, or restore the integrity, connectivity, and resilience by increasing the area of natural ecosystems by 2050. Regarding forest ecosystems, the extent area is a headline indicator, combined with complementary measures such as forest distribution, tree cover loss, forest fragmentation index, forest landscape integrity index, and biomass density (CBD-UN, 2022). The SEEA-EA suggests a wide range of biophysical indicators to determine forest ecosystem conditions (United Nations et al., 2021). For instance, the structural status is defined by Aboveground biomass (AGB) density and tree cover density (TCD); the functional status can be represented by net primary productivity (NPP) and dry matter productivity (DMP); the compositional state indicators capture ecosystem diversity, species and abundance; and connectivity and fragmentation indicators measured through landscape characteristics, such as number and size of patches, and other several shape and edge indices. AGB, TCD, NPP, and DMP are assessed via ground measurements, although these are expensive and challenging to deploy, especially across large areas (Paul et al., 2022). Remote sensing data and technologies strengthen ecosystem condition mapping, providing continuous observation at a large scale with relatively low costs (Rapinel et al., 2015). The advancement of Geographical Information Systems (GIS) and power computation has provided a wealth of analytical tools for comprehensively analyzing landscape characteristics and the relationship between the environment and society (Mengist et al., 2021).

Landscape metrics offer a means to depict the composition and organization of a landscape, providing insights into the arrangement and diversity of its elements (Haines-Young et al., 2012). These quantitative measures capture the landscape's spatial structure, distribution patterns, dynamics, and overall evolution (Holzwarth et al., 2020). Recent improvements and innovations in predictive land-use models have led to significant advances in spatial modelling, capturing complexities of land-use systems, analyzing drivers of changes and patterns, and simulating future scenarios (Roodposhti et al., 2020). Predicting future dynamics using historical land use data has proven to be a valuable mechanism, not only for anticipating future trends but also for enhancing and deepening the understanding of past dynamics (Anselmetto et al., 2022). Dynamic modelling based on cellular automata (CA) is one of the methods developed for forecasting environmental patterns and changes (Moreno et al., 2010) providing potential outcomes for land-use policy scenarios (de Brito et al., 2021). CA models functioning basis on a discrete dynamic system that divides space into regular spatial cells and moves time in discrete increments according to predefined rules (Hewitt et al., 2019).

Forests have been widely studied on a global scale, yet there remains a need for advanced methods able to provide consistent data to better understand the interactions between forest resources and human activities (Cunha et al., 2021; Holzwarth et al., 2020; Zarandian et al., 2018). For instance, research has focused on various aspects of forest dynamics, such as estimating the risk of forest-to-non-forest conversions (Fitts et al., 2021; Thanh Noi and Kappas, 2017) and vulnerabilities due to fire, logging, climate change, and pests (Chave, 1999; Franklin et al., 2002; Houghton et al., 2009). In addition to building knowledge regarding the consequences of urban sprawl on forest ecosystems and biodiversity (Mengist et al., 2021; Mulatu et al., 2017; Paul et al., 2022), the fragmentation effects due to changes in the structure of the ecological landscape (Bossel and Krieger, 1991; Dou et al., 2023; Hoffmann et al., 2022; Madrigal-González et al., 2023) and the negative impacts on natural capital and ES (Golub et al., 2009; Holzwarth et al., 2020). Despite, the numerous studies mapping forest dynamics, there is still insufficient research about present and future forest fragmentation impacts, and acknowledging the most relevant indicators of forest condition. Thus, the objective of this study is to provide a list of relevant measurements that can be systematically acquired, combining satellite-derived data and computational modelling, serving as complement indicators to decision support frameworks such as the SEEA-EA and the CBD-UN. In this scope, the article focuses on forest ecosystems including coniferous, broad-leaved forests and mixed forests, occurring at the country level of Mainland Portugal. The methodology includes reclassifying the Corine Land Cover maps of years 1990, 2000, 2006, 2012, and 2018, used as reference data to predict national forest distribution for five decades utilizing a CA algorithm through the SIMLANDER software (Hewitt et al., 2019). A total of 130 landscape metrics (LM) is calculated using the landscape metrics package (Hesselbarth et al., 2019), filtered out through a variance threshold method and the top 10 indicators were selected and ranked using dimensionality reduction and feature importance techniques from scikit-learn modules (Pedregosa et al., 2011). The overall outcome of this study will deepen the knowledge of the spatial influence of forest configuration and how future scenarios could be modelled considering anthropogenic and environmental pressures. Additionally, recognizing the most effective LM-based indicators of forest condition helps to more accurately monitor ES, plan
landscape development and resource management, prevent impacts on
biodiversity and natural capital, and track progress towards sustainable
development.

**2. Materials and methods**

**2.1. Study area**

Portugal is set in the Mediterranean Basin, one of the five regions of the world characterized by Mediterranean-type ecosystems. This biogeographic region is known for its climatic regimes of mild and wet winters and warm and dry summers (Marques et al., 2019). The annual precipitation averages roughly 3000 mm in the north and 500 mm in the south, with substantial seasonality (de Lima et al., 2015). Spring (March through May) accounts for 24 % of annual precipitation, summer (June through August) accounts for 6 %, fall (September through November) accounts for 28 %, and winter (December through February) accounts for 42 % (Coelho et al., 2013).

The low precipitation occurring in the summer favours forest fires as vegetation becomes quickly dry and flammable but also increases drought severity, deforestation, rural abandonment, and water scarcity (Chrysafis et al., 2020). Seasonal fluctuations alter the production of vegetation that depends on the presence of freshwater or subterranean water (Wang et al., 2023). Due to restricted water and nutrient availability, it often peaks during the wet season and declines during the dry season (Mpakairi et al., 2022). The combination of higher latitude and altitude and the annual rainfall pattern in the Northern region create better environmental conditions for forest ecosystems than in the Southern region (Belo-Pereira et al., 2011). Due to low rainfall, the south of Portugal is known as a water-scarce region potentially exposed to climate change and its environmental, social, and economic consequences (Almeida and Cabral, 2021). Forests near lakes, springs, and perennial streams take advantage of water availability during the reproductive period, which matches with the hottest season of the year (Harris, 1984), and large trees can store significant amounts of intra- and intercellular water during the winter which can be used for growth the following summer (Wilcox and Murphy, 1985). Forest and semi-natural areas and croplands are the most represented land use classes at the

national level (Copernicus Programme, 2023). The flat terrain of the southern regions is dominated by extensive agriculture and pastures integrated into agro-forestry systems, where low-density evergreen oak trees like Cork oak and Holm oak are prevalent (Santos et al., 2023). In contrast, the higher elevations in the central and northern regions are marked by small landholdings, family-run farms, and forests primarily composed of pine and eucalyptus (Fonseca et al., 2019). Fig. 1A) displays an overview of the study area within Europe and B) the Land Use and Land Cover Map of the year 2018 (Copernicus Programme, 2023).

**2.2. Methods**

The methodology was structured in three main parts: 1) predicting forest dynamics and evaluating models' accuracy through quantity and allocation disagreements; 2) calculating 130 LM, filtering them out, selecting and ranking the top 10 indicators; and 3) assessing forest ecosystem conditions for Mainland Portugal. Fig. 2 summarizes the methodological steps, specifying the software utilized, implemented tasks, input datasets and outcomes.

**2.2.1. Modelling forest dynamics**

The process starts by mapping the presence/absence of forests at the national level by reclassifying the CORINE land cover maps of 1990, 2000, 2006, 2012, and 2018. Land use and land cover (LULC) datasets were obtained from the Copernicus Land Monitoring Services (Copernicus Programme, 2023). Each raster was reclassified as a binary map of forest/non-forest in ArcGIS Pro version 3.3 (ESRI, 2023). The forest class comprises broad-leaved, coniferous, and mixed forests, with a vegetation pattern composed of native or alien coniferous and/or broad-leaved trees capable of producing timber or other forest products. The forest trees should be able to reach a minimum height of 5 m, canopy closure is greater than 30 %, and the minimum threshold for young plantations is 500 trees ha⁻¹ (European Environment Agency, 2023).

Spatial predictions can be modelled through open-source software such as the SIMLANDER (Hewitt et al., 2022), a constrained cellular automata (CA) system that simulates future scenarios of land use based on land claims (Hewitt et al., 2019). CA models are discrete dynamical

<figure><img src="image_6.png" /><figcaption>**Fig. 1.** Study area: A. Location of Mainland Portugal in Europe; B. Corine Land Use and Land Cover Map – Portugal 2018 (Copernicus Programme, 2023).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 2.** Summary of methodological steps, specifying the software utilized, implemented tasks, input data and outcomes. The methodology is organized into three main tasks: 1) predict forest dynamics; 2) select the top 10 indicators; and 3) assess forest ecosystem conditions for Mainland Portugal.</figcaption></figure>

systems (in space and time) whose behaviour is specified in terms of a local relation defined by rules (Jokar Arsanjani et al., 2013). The transition rules are fundamental in such models as they are responsible for the state of the cell, at each time step; they quantify the spatial effect that the forecast cells contain in the LULC changes (Ozturk, 2015). SIMLANDER version 2.0.0 can take two land use maps on different dates as input with the rationale that the future cell state will be defined by the

state of the same cell and immediate neighbouring cells at the present cell state (Moreno et al., 2010). Table 1 lists datasets, respective sources, and specifications of use and application for the SIMLANDER software.

The simulation approach also considers other drivers known to influence LULC change, such as the degree of attraction of a particular land use to others in its neighbourhood (N) and the distance at which they are attractive, the proximity to infrastructure networks, known as accessibility (A), the underlying land suitability for development (S) and a random factor (v) to account for unknown or undetermined factors that influence real-world outcomes (Hewitt et al., 2020). The national infrastructure networks were obtained from the OpenStreetMap project (OpenStreetMap, 2023), and the topography from the National Aeronautics and Space Administration (NASA) (NASA, 2023) through the Shuttle Radar Topography Mission (SRTM). The elevation dataset was used to calculate the slope, and the infrastructure networks were used to calculate the Euclidian distance of forests to highways and railways. The accessibility map referred to the distance (Euclidian) of forests to highways and railways, and the suitability map associated with the capability of forests to be developed in specific locations based on the characteristics of the terrain, were processed and resampled through interpolation (nearest neighbour assignment) using the Resample tool of the ArcGIS Pro. Both rasters were upscaled to the same spatial resolution of the reference land-use maps (100 m).

The number of runs was set to 50, and the option for a new random seed factor was selected to be generated each time. The number of cells allocated in the model is controlled by the cell demand, which is determined by dividing the number of new cells of non-forest land that have emerged between the start and end dates and the number of years between them (Hewitt et al., 2022). We have estimated the demand values for scenarios 2036 and 2054 through linear interpolation, and the

<caption>**Table 1**</caption>

Description of the datasets used as input in the SIMLANDER software version 2.0.0 (Hewitt et al., 2022), including their source, parametrization, and application.

<table><thead><tr><td>Dataset</td><td>Source</td><td>Model parameter</td><td>Application</td></tr></thead><tbody><tr><td>CORINE Land use land cover maps (100 m)</td><td>Copernicus Land Monitoring Services (https://land.copernicus.eu/)</td><td>Initial and final maps</td><td>Used as reference maps and for the assessment of models' performance</td></tr><tr><td>Topography (30 m)</td><td>National Aeronautics and Space Administration (NASA) - Shuttle Radar Topography Mission (https://www.earthdata.nasa.gov/)</td><td>Suitability map</td><td>Used to characterize morphological aspects of the terrain</td></tr><tr><td>Infrastructure networks (50 m)</td><td>OpenStreetMap project (https://www.openstreetmap.org/)</td><td>Accessibility map</td><td>Used to calculate the Euclidian distance within forests to highways and railways</td></tr></tbody></table>
models' parameters were tested and analysed to calibrate the predictions and select the best models. Table S1 in the Supplementary Materials lists simulations by year, describing input data and the model's parameters for the best-selected model.

Measuring prediction quality and accuracy is one of the most critical steps in geospatial modelling (Mengist et al., 2021). The accuracy and reliability of predictions from CA models depend on model design, appropriate parameterization, validation against empirical data, and consideration of uncertainties, due to data limitations, model assumptions, and spatial and temporal resolution (Roodposhti et al., 2020). The CA predictions were assessed by accuracy functions including overall accuracy (OA), user's accuracy (UA), producer's accuracy (PA), recall, precision, and F1-score. The accuracy assessment workflow in ArcGIS Pro version 3.3 respectively employs the following tools: Create Accuracy Assessment Points, Update Accuracy Assessment Points, and Compute Confusion Matrix. The first tool creates a set of random points and assigns a class to them based on the reference data (reclassified Corine LULC maps - forest/non-forest classes). The second tool updates the target field in the attribute table to compare reference points to the CA prediction maps. In both tools, the total number of random points generated equals 500,000, and the sampling strategy was the Stratified Random option (points randomly distributed within each class proportional to its relative area). Both tools ensure that each point will have valid class values for the predictions and reference fields enabling the third tool to compute the confusion matrix. The recall and precision were further manually calculated based on the results of UA, PA, OA, and F1-scores. This workflow had as input the observed and predicted maps of the years 2000, 2006, 2012, and 2018.

The UA function is related to false positives (FP), or errors of commission, where pixels incorrectly classified as forest are not forest. UA is calculated by dividing the number of correctly predicted forest points by the total number of points predicted as forest. PA deals with false negatives (FN), or errors of omission, indicating how well the prediction results align with the reference data, and it is calculated by dividing the number of correctly predicted points by the total number of reference points for that class. OA represents the proportion of reference points that were correctly classified out of the total number of reference points. Recall is the ratio of true positives (TP) to the sum of TP and FN. It reflects the model's ability to identify all forest samples. Precision is the ratio of TP to the sum of TP and FP. It measures the model's ability to avoid incorrectly labelling non-forest samples as forest. The F1-score is the harmonic mean of precision and recall, providing a balanced measure of the model's accuracy. All these functions follow the convention of higher values and better predictions.

While these functions provide useful information on model quantity disagreements, they do not explicitly address the magnitude to which simulations match the spatial arrangement of observed landscape, such as patterns of fragmentation, clustering, or connectivity within patches (Pickard and Meentemeyer, 2019). Analyzing the spatial pattern matching through assessing the betweenness, degree and closeness of forest patches goes beyond quantifying disagreements (Foody, 2020), as they measure the importance, influence, and accessibility of a node (the centroid of each forest patch) within a network (landscape) (Gabriel, 1969). The betweenness measures how often a node lies on the shortest path between other nodes in the network (Matula and Sokal, 2010). The range of values depends on the structure of the network, the minimum value is 0, occurring when a node does not lie on any shortest paths between other nodes (Savary et al., 2021). The maximum possible value depends on the network size and topology, and it occurs when a node lies on all shortest paths between every pair of other nodes (Matula and Sokal, 1980). The degree tells us about the immediate connectivity of a patch and how many direct connections it has (Csardi and Nepusz, 2006). It gives a quick sense of the local importance of a node, as it reflects how well-connected a node is within its immediate network (Freeman, 1978). The minimum value is 0, occurring when a node is isolated, meaning it has no connections to any other nodes; and the

maximum possible value occurs when a node is connected to all other nodes in the network (Marchette, 2004). Closeness measures how close a node is to all other nodes in the network, in terms of the shortest path distance (Gabriel, 1969). The smallest value occurs when a node is far from all other nodes. For an unconnected or highly distant node, the closeness approaches 0; and the maximum value is 1, which occurs when the node is directly connected to every other node (Brandes, 2001). The landscape graphs construction and analyses were conducted through the graph4lg v. 1.8.0 package in R (Savary et al., 2021).

The assessment of forest ecosystem conditions followed the SEEA-EA guidelines on Biophysical Modelling for Ecosystem Accounting (Edens et al., 2022), which includes the AGB of forests as an indicator of ecosystem condition, and a proxy for productivity, and carbon stock and sequestration (King et al., 2024). The mean AGB density (110.50 tons ha⁻¹) value for mainland Portugal was provided by the Institute for Nature Conservation and Forests, I.P. (ICNF, 2023) and converted into carbon stock assuming that 50 % of the dry AGB of forest corresponds to carbon (Chave, 1999; Goetz et al., 2009; Hojo et al., 2023). To predict the national carbon stock for the studied years we multiplied these values by the estimated forest extents.

### 2.2.2. Landscape analysis

Quantitative grouping of similar landscape change trends is an important part of landscape ecology due to the relationship between a pattern and an underlying ecological process (Nowosad and Stepinski, 2019a). Spatial and temporal variations in LULC capture different components of landscapes that are considered important for assessing biodiversity and ES (Fletcher and Fortin, 2018), describing landscapes at the patch, class, and landscape levels. A patch is defined as neighbouring cells belonging to the same LULC class. Class-level metrics describe all patches belonging to the same class, and lastly, landscape-level metrics describe the whole landscape composed of all patches (McGarigal et al., 2012). For analysis, we used the landscapemetrics package v.2.1.1 (Hesselbarth et al., 2019) developed through the programming language R 3.6.0 (R Core Team, 2021). The package reimplments the most common metrics from FRAGSTATS 4.0 (McGarigal and Marks, 1995) and new ones that have been further developed, such as the marginal entropy, conditional entropy, joint entropy, and mutual information included in the information-theoretical framework. The metrics are grouped into six different types according to the characteristic of the landscape they describe: 1) area and edge metrics describe the size of patches and classes and the amount of edge, characterizing the composition of the landscape and dominance of classes; 2) shape metrics describe the shape of patches, mainly by using its area and perimeter; 3) core metrics refers to the portion of a patch that is not part of the edge, denoting the area within a patch where the pixels are sufficiently distant from the boundary; 4) aggregation metrics not only describe the degree of aggregation or separation of patches but also capture patterns of interspersion and intermixing, defining the distribution and spatial arrangement of patches in relation to their neighbouring patches; 5) diversity metrics accounts for both the number of patch types and their relative abundance, and how the patches are distributed across different land-cover types; and 6) complexity metrics are applied to quantification and classification of landscape patterns, reflecting the degree of similarity or correlation between different trends. Since there is no consensus in the literature on which metrics are the most appropriate for this specific problem (Roodposhti et al., 2020), and including those highlighted by SEEA-EA, such as the number of patches and mean patch size, Simpson's diversity index, shape and edge areas, all metrics available on the package (a total of 130) were calculated. The metrics are grouped at the landscape level (a total of 64); at the class level (54); and the patch level (12). Few metrics exist within all levels, except the diversity metrics that are only available on the landscape level. The only metric that could not be determined was the iji (Interspersion and Juxtaposition index - Aggregation metric), which requires a number of classes larger than three.
The selection process started calculating summary statistics and removing features whose variance does not meet a certain threshold using the Variance Threshold algorithm, an unsupervised feature selection module from the scikit-learn v.1.5.1 (Pedregosa et al., 2011). The threshold setting was tested by analyzing the results of forest variation percentage (1) to remove a maximum number of features without losing information.

$$ \text{Forest variation (\%)} = \frac{\text{metric}_{2018} - \text{metric}_{2000}}{\text{metric}_{2000}} \times 100 \quad (1) $$

The following method was the Principal Component Analysis (PCA) from the decomposition module, combined with an unsupervised Feature Importance algorithm from the feature selection module. The PCA function looks for a combination of attributes (principal components, or directions in the feature space) that account for the most variance in the data (Maćkiewicz and Ratajczak, 1993). It reduces the number of features while preserving as much of the original information as possible, transforming the original variables into a new set of uncorrelated variables called principal components (Elhaik, 2022). After calculating the loadings for the first principal component (PC1), the features were sorted by their relevance using the Select from Model algorithm. The scripts used for the analysis and selection of LM-based indicators are detailed in Table S2 (Supplementary Materials).

# 3. Results

## 3.1. Evaluation of forest dynamic predictions

Based on past trends from 1990 to 2018 (Fig. 3A), the predictions estimated continued forest loss, even under the most optimistic scenarios (Fig. 3B). In the short-term projection (2036), the simulation shows moderate changes in forest cover, such as gradual forest loss, expansion of non-forest areas, and increased forest fragmentation. For 2054 the map exhibits more pronounced changes, including larger areas of deforestation and shifts in forest distribution.

The accuracy assessment workflow evaluated OA, PA, UA, recall, precision, and F1-score (Table 2) for forest maps across four years: 2000, 2006, 2012, and 2018. All accuracy functions show a significant drop from 2000 to 2018. While the decline from 2000 to 2006 is sharp, the metrics stabilised somewhat after 2006, suggesting that the forest predictions did not degrade as drastically in the later years. The OA was highest in 2000 (96.6 %) but dropped significantly in 2006 (87.8 %) and stabilizes in 2012 and 2018 at around 86–87 %. The decline in PA is stark, going from 97.7 % in 2000 to 64 % in 2018. A similar trend is observed for UA, recall and precision, with a drop from 93 % in 2000 to around 64 % in 2018. The F1-score, which balances precision and recall, follows the same trend as the other metrics. It started high in 2000 at 93.6 % and dropped sharply to around 64 % by 2018.

The number of forest patches on the reference maps increased 19.3 % from 2000 to 2018 (Table 3). For predictions, it shows a declining trend over time of 42.5 % from 2000 to 2054. The observed mean betweenness value increases over time from 1.5*10⁵ in 2000 to 2*10⁵ in 2018, in contrast, predicted values decrease from 2*10⁵ in 2000 to 9.3*10⁴ in 2054. The observed mean degree remains relatively stable over time, fluctuating slightly between 7.8 and 8, and for predictions values range from 7.7 to 7.9. The observed mean closeness values decreased from 7.1*10⁻² in 2000 to 6.5*10⁻² in 2018. The predicted values show the

<table><caption>**Table 2**
Performance evaluation of the SIMLANDER model for forest dynamic predictions. The accuracy assessment workflow was conducted in ArcGIS Pro v.3.3: Overall accuracy (OA), producer's accuracy (PA), user's accuracy (UA), recall, precision, and F1-score.</caption><thead><tr><th>Forest Map (year)</th><th>OA (%)</th><th>PA (%)</th><th>UA (%)</th><th>Recall (%)</th><th>Precision (%)</th><th>F1-score (%)</th></tr></thead><tbody><tr><td>2000</td><td>96.6</td><td>97.7</td><td>93.7</td><td>93.7</td><td>93.6</td><td>93.6</td></tr><tr><td>2006</td><td>87.8</td><td>72.4</td><td>72.4</td><td>72.4</td><td>72.4</td><td>72.4</td></tr><tr><td>2012</td><td>86.4</td><td>67.4</td><td>67.4</td><td>67.4</td><td>67.4</td><td>67.4</td></tr><tr><td>2018</td><td>86.7</td><td>64</td><td>64.2</td><td>64.2</td><td>64</td><td>64.1</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Fig. 3.** Forest distribution in mainland Portugal. A. Reference maps used as input for modelling future projections of forest/non-forest areas; B. Predicted maps of forest obtained from the SIMLANDER model.</figcaption></figure>
<caption>**Table 3**</caption>

Number of forest patches at the national scale, and values of mean betweenness, mean degree, and mean closeness calculated in each patch centroid. These values were obtained from the graph4lg v. 1.8.0 package in R (Savary et al., 2021). Legend: obs - observed, pred - predicted. The future projections column highlights estimations for 2036 and 2054.

<table><thead><tr><th>Indicators</th><th>2000</th><th>2006</th><th>2012</th><th>2018</th><th>2036</th><th>2054</th><th>Future projections</th></tr></thead><tbody><tr><td>Number of patches obs</td><td>5.4*10<sup>3</sup></td><td>6.3*10<sup>3</sup></td><td>6.4*10<sup>3</sup></td><td>6.7*10<sup>3</sup></td><td></td><td></td><td></td></tr><tr><td>Number of patches pred</td><td>6.6*10<sup>3</sup></td><td>6.5*10<sup>3</sup></td><td>4.7*10<sup>3</sup></td><td>4.3*10<sup>3</sup></td><td>3.6*10<sup>3</sup></td><td>3.8*10<sup>3</sup></td><td>x</td></tr><tr><td>Mean betweenness obs</td><td>1.5*10<sup>5</sup></td><td>1.9*10<sup>5</sup></td><td>1.8*10<sup>5</sup></td><td>2*10<sup>5</sup></td><td></td><td></td><td></td></tr><tr><td>Mean betweenness pred</td><td>2*10<sup>5</sup></td><td>1.9*10<sup>5</sup></td><td>1.2*10<sup>5</sup></td><td>1.1*10<sup>5</sup></td><td>8.3*10<sup>4</sup></td><td>9.3*10<sup>4</sup></td><td>x</td></tr><tr><td>Mean degree obs</td><td>7.8</td><td>8</td><td>8</td><td>7.9</td><td></td><td></td><td></td></tr><tr><td>Mean degree pred</td><td>7.7</td><td>8</td><td>7.9</td><td>7.8</td><td>7.8</td><td>7.8</td><td>x</td></tr><tr><td>Mean closeness obs</td><td>7.1*10<sup>-2</sup></td><td>6.7*10<sup>-2</sup></td><td>6.6*10<sup>-2</sup></td><td>6.5*10<sup>-2</sup></td><td></td><td></td><td></td></tr><tr><td>Mean closeness pred</td><td>6.3*10<sup>-2</sup></td><td>6.6*10<sup>-2</sup></td><td>7.7*10<sup>-2</sup></td><td>8.1*10<sup>-2</sup></td><td>8.6*10<sup>-2</sup></td><td>8.6*10<sup>-2</sup></td><td>x</td></tr></tbody></table>

opposite trend, with mean closeness increasing from $6.3*10^{-2}$ in 2000 to $8.6*10^{-2}$ in 2054. As the future simulations were optimistically calibrated, the results showed slightly opposite value patterns. The spatial pattern analysis based on forest patch centroids was conducted through Gabriel's graphs (Gabriel, 1969) (Fig. 4 - betweenness (A), degree (B), and closeness (C)).

Following the SEEA-EA framework, we quantified the national physical ecosystem accounts, by assessing forest extent and indicators of

its condition considering observations and predictions (Table 4). The results evidence decreased forest area between 2000 and 2018, representing forest area total losses of 30.5 % equivalent to $7.3*10^5$ ha. Consequently, the national AGB carbon stock decreased from 131.5 to 91.3 Megatons (Mt), i.e. a decrease of 31 % for the observed maps, against a reduction of 32.8 % for predictions.

<figure><img src="image_9.png" /><figcaption>**Fig. 4.** Gabriel graphs for the spatial pattern analysis of betweenness (A), degree (B), and closeness (C) in each study year. In (A) the range of values depends on the structure of the network, and the minimum value is 0, when a node (forest patch centroid) does not lie on any shortest paths between other nodes, and is maximum when a node lies on all shortest paths between every pair of other nodes. For (B) the minimum possible value is 0, occurring when a node is isolated, meaning it has no connections to any other nodes; and the maximum possible value occurs when a node is connected to all other nodes in the network. In (C) the smallest value occurs when a node is far from all other nodes. For an unconnected or highly distant node, the closeness approaches 0; and the maximum value is 1, which occurs when the node is directly connected to every other node. Legend: Obs - observed; Pred - predicted.</figcaption></figure>
<caption>**Table 4**</caption>

National forest ecosystem extent and AGB carbon (AGB C) stock for observed and predicted maps. The mean forest AGB density at the national scale (110.50 tons ha⁻¹) was provided by the Institute for Nature Conservation and Forests of Portugal (ICNF, 2023). This value was used as a proxy to estimate carbon stock for the studied years, based on ecosystem extent and considering that 50 % of the dry AGB of forest corresponds to carbon (Chave, 1999; Goetz et al., 2009; Hojo et al., 2023). Legend: obs - observed, pred - predicted. Future projections refer to the calculations that were predicted for 2036 and 2054.

<table><thead><tr><th>Indicators</th><th>2000</th><th>2006</th><th>2012</th><th>2018</th><th>2036</th><th>2054</th><th>Future projections</th></tr></thead><tbody><tr><td>Extent obs (ha)</td><td>2.3*10<sup>6</sup></td><td>1.3*10<sup>6</sup></td><td>1.5*10<sup>6</sup></td><td>1.6*10<sup>6</sup></td><td></td><td></td><td></td></tr><tr><td>Extent pred (ha)</td><td>2.3*10<sup>6</sup></td><td>1.4*10<sup>6</sup></td><td>1.5*10<sup>6</sup></td><td>1.6*10<sup>6</sup></td><td>1.3*10<sup>6</sup></td><td>1*10<sup>6</sup></td><td>x</td></tr><tr><td>AGB C stock obs (Mt)</td><td>131.4</td><td>75.5</td><td>83</td><td>91.3</td><td></td><td></td><td></td></tr><tr><td>AGB C stock pred (Mt)</td><td>131.5</td><td>75.3</td><td>82.8</td><td>88.4</td><td>71.8</td><td>55.3</td><td>x</td></tr></tbody></table>

## 3.2. Landscape metrics analysis

The first step of filtering out features includes testing different thresholds of variance, and identifying evolutional trends of forest dynamics calculating the mean percentage variation. The variance threshold column lists the values used as thresholds in the tests, for instance, for threshold values between 0.05 and 0.3, and 0.4 to 0.9, the number of features removed and the mean variation percentage within these intervals were the same. Considering the variance threshold that removed the maximum number of features (values between 0.4 and 0.9) we selected the features that were above the highest mean variation (50 %) within this interval, corresponding to 27 features (Table 5). These metrics are presented in detail and grouped by the type of metrics in the Supplementary Materials (Table S3).

The selected 27 LM were then analysed through PCA by 10 components (Fig.S1 – Supplementary Materials). The explained variance percentages for PC1, PC2 and PC3 are respectively 62.5 %, 18.4 % and 16.2 %. The top 10 loadings for PC1, are the features that contribute the most to this principal component and capture the largest variance in the dataset. These 10 most important features were selected as LM-based indicators of forest condition (Table 6).

The effective mesh size contributes the most to PC1, suggesting a strong influence on the main source of variation in the data. It measures the degree of landscape fragmentation by considering patch size distribution. A favourable change would be an increase in its measure indicating lower fragmentation and higher landscape connectivity. The second most important metric is the splitting index indicating that areas with more fragmentation (splitting) heavily contribute to the variation captured by PC1. A favourable change would be a decrease in the value, indicating less fragmentation and larger patches. The mean radius of gyration is the third most representative metric describing the extent to which patch areas are spread out from their centroids, reflecting the landscape connectivity. A favourable change is an increase in the measurements indicating that patches are more spatially connected. Next, is the largest patch index describing the percentage of the total landscape area occupied by the largest patch. A favourable change is an increase in its measurement suggesting a dominant patch is getting larger. The fifth metric most relevant to PC1 is the mean of the core area showing

favourable changes when its values increase as a result of reduced fragmentation and increased patch sizes. A positive loading suggests that larger core areas (unfragmented central parts of patches) are heavily correlated with PC1, indicating a relationship between core size and overall landscape variance. Next is the core area percentage of landscape, which a favourable change is an increase in the measure indicating a rise in the availability of forest habitat in the landscape. A high loading means that PC1 captures variation related to the total area covered by forest patches, indicating that areas with more total patch area are significant contributors to the variance. The Simpson's evenness index is the seventh metric contributing the most to PC1, its positive contribution means that more fragmented or isolated patches are strongly represented in PC1, evidencing that PC1 is capturing fragmentation patterns. The mutual information reflects the degree of similarity or correlation between different landscape patterns. A favourable change is an increase in its value suggesting greater similarity between patterns, implying better structural connectivity and coherence in landscape configuration. A strong loading indicates that PC1 is capturing landscape complexity. The Simpson's diversity index accounts for both the number of patch types and their relative abundance. An increase in this index is considered a favourable change, and a high positive loading suggests that areas with greater landscape diversity are significant contributors to the variance captured by PC1. The tenth select metric is the mean of the contiguity index representing patch contiguity (how connected patches are) and cohesion within individual patches. A high loading means that more contiguous patches, or landscapes with higher connectivity, are heavily represented in PC1. A favourable change is an increase in the value meaning that patches are more compact or contiguous.

These selected indicators were then used to assess the forest ecosystem conditions in Mainland Portugal, for both normalized observations (Fig. 5A) and predictions (Fig. 5B). These results are presented in detail in the Supplementary Materials (Table S4).

Regarding the observations the mesh value decreases 13 % between 2000 and 2018, indicating a significant increase in landscape fragmentation. For predictions, the results follow a similar pattern. Within observed maps, the mean core area decreased by 48.7 % from 2000 to 2018, while for the predictions it did not fluctuate that much a part for the year 2012. The mean radius of gyration decreased within observations but increased within predictions. The split index for the reference maps increased about six times more from 2000 to 2018, and for predictions, the index shows a similar rising trend. The core area percentage decreased from 18.9 % in 2000 to 11.9 % in 2018 in the observed maps, while the predictions show a steady similar decline, from 19.5 % in 2000 to 8.9 % in 2054. The largest patch index decreased from 5.7 % in 2000 to 1.4 % in 2018, according to the observations, and a similar decline in predictions, dropping from 5.6 % in 2000 to 0.9 % in 2054. The shape index decreased from 0.8 in 2000 to 0.6 in 2018 in the observations, and from 0.8 in 2000 to 0.4 in 2054 for the predictive maps. The mean contiguity index remains stable at 0.7 throughout 2000–2018 for the observations, and oscillates between 0.6 and 0.7 for the predictions, indicating that, while fragmentation increases, the patches that remain tend to be relatively compact and contiguous. The mutual information index remains stable at 0.5 from 2000 to 2018 in the observations but

<caption>**Table 5**</caption>

Variance threshold tests to filter out features with zero or low variances. The variance threshold column lists the values used as thresholds in the tests, and the mean variation percentage column depicts the value within the number of features remaining in some threshold.

<table><thead><tr><th>Variance threshold</th><th>Number of features removed</th><th>Number of features remained</th><th>Mean Variation %</th><th>Number of features remained within the variance threshold and above the mean variation</th></tr></thead><tbody><tr><td>0.01</td><td>45</td><td>85</td><td>13 %</td><td>66</td></tr><tr><td>0.05; 0.1; 0.2; 0.3</td><td>51</td><td>79</td><td>36 %</td><td>49</td></tr><tr><td>0.4; 0.5; 0.6; 0.9</td><td>53</td><td>77</td><td>50.3 %</td><td>27</td></tr></tbody></table>
<table><caption>**Table 6**</caption><thead><tr><th>Name (unit)</th><th>Type</th><th>Level</th><th>Formula, Range [min, max], Reference</th><th>Favourable change</th><th>Importance to PC1</th></tr></thead><tbody><tr><td>Effective mesh size (ha)</td><td>Aggregation</td><td>class</td><td>mesh = <sup>&sum;<sub>j=1</sub><sup>n</sup> a<sub>ij</sub><sup>2</sup></sup>&frasl;<sub>A</sub> &times; <sup>1</sup>&frasl;<sub>10000</sub><br>A is the total landscape area, a<sub>ij</sub> is the patch area<br>Range [cell size/A, A]; mesh equals min value whenever the class covers only one cell and equals max if only one patch is present.<br>(Jaeger, 2000; McGarigal et al., 2012)</td><td>Increase</td><td>2.5*10<sup>-1</sup></td></tr><tr><td>Splitting index (none)</td><td>Aggregation</td><td>class</td><td>split = <sup>A<sup>2</sup></sup>&frasl;<sub>&sum;<sub>j=1</sub><sup>n</sup> a<sub>ij</sub><sup>2</sup></sub><br>A is the total landscape area, a<sub>ij</sub> is the patch area<br>Range [1, number of cells squared]; split equals min value if only one patch is present. It reaches the max as the number of patches increases.<br>(Jaeger, 2000; McGarigal et al., 2012)</td><td>Decrease</td><td>2.4*10<sup>-1</sup></td></tr><tr><td>Mean radius of gyration (m)</td><td>Area and edge</td><td>landscape</td><td>gyrate_mn = mean(<sup>&sum;<sup>z</sup><sub>r=1</sub></sup>&frasl;<sub>z</sub> h<sub>ijr</sub> [patch<sub>ij</sub>])<br>h<sub>ijr</sub> is the distance from each cell to the centroid of each patch and z is the number of cells<br>Range [0, inf[; gyrate_mn equals min whenever a patch is a unique cell. It reaches the max range when only one patch is present.<br>(Keitt et al., 1997; McGarigal et al., 2012)</td><td>Increase</td><td>2.4*10<sup>-1</sup></td></tr><tr><td>Largest patch index (%)</td><td>Area and edge</td><td>class</td><td>lpi = <sup>n</sup>&frasl;<sub>j=1</sub> max(a<sub>ij</sub>) &times; <sup>100</sup>&frasl;<sub>A</sub><br>A is the total landscape area, a<sub>ij</sub> is the patch area<br>Range [0,100]; lpi equals min when the largest patch is becoming small and equals max when only one patch is present.<br>(McGarigal et al., 2012)</td><td>Increase</td><td>2.4*10<sup>-1</sup></td></tr><tr><td>Mean of core area (ha)</td><td>Core area</td><td>class</td><td>core_mn = mean(a<sub>ij</sub><sup>core</sup> [patch<sub>ij</sub>])<br>a<sub>ij</sub><sup>core</sup> is the core area<br>Range [0, inf[; core_mn equals min if core equals 0 for all patches. It reaches max as the core area increases.<br>(McGarigal et al., 2012)</td><td>Increase</td><td>2.4*10<sup>-1</sup></td></tr><tr><td>Core area percentage of landscape (%)</td><td>Core area</td><td>class</td><td>cpland = (<sup>&sum;<sup>n</sup><sub>j=1</sub></sup>&frasl;<sub>A</sub> a<sub>ij</sub><sup>core</sup>) &times; 100<br>a<sub>ij</sub><sup>core</sup> is the core area and A is the total landscape area<br>Range [0, 100]; core equals 0 for all patches. It reaches max as the amount of core area increases.<br>(McGarigal et al., 2012)</td><td>Increase</td><td>2.4*10<sup>-1</sup></td></tr><tr><td>Simpson's evenness index (none)</td><td>Diversity</td><td>landscape</td><td>siei = <sup>1 - &sum;<sub>i=1</sub><sup>m</sup> P<sub>i</sub><sup>2</sup></sup>&frasl;<sub>1 - <sup>1</sup>&frasl;<sub>m</sub></sub><br>P<sub>i</sub> is the proportion of class i and m is the number of classes<br>Range [0,1]; siiei is min when only one patch is present. It reaches max when the number of classes increases within equally distributed proportions.<br>(May, 1975; McGarigal et al., 2012; Romme, 1982; Simpson, 1949)</td><td>Increase</td><td>2.4*10<sup>-1</sup></td></tr><tr><td>Mutual information (none)</td><td>Complexity</td><td>landscape</td><td>mutinf = I(y,x)<br>I is the information, y and x are random variables<br>Range [0,1[; mutinf increases with landscape diversity.<br>(Nowosad and Stepinski, 2019b)</td><td>Increase</td><td>2.4*10<sup>-1</sup></td></tr><tr><td>Simpson's diversity index (none)</td><td>Diversity</td><td>landscape</td><td>sidi = 1 - &sum;<sub>i=1</sub><sup>m</sup> P<sub>i</sub><sup>2</sup><br>P<sub>i</sub> is the proportion of class i<br>Range [0,1[; sidi is min when only one patch is present. The value increases as the number of class types increases within equally distributed proportions.<br>(May, 1975; McGarigal et al., 2012; Romme, 1982; Simpson, 1949)</td><td>Increase</td><td>2.4*10<sup>-1</sup></td></tr><tr><td>Mean of contiguity index (none)</td><td>Shape</td><td>landscape</td><td>contig_mn = mean(<sup>(<sup>&sum;<sup>z</sup><sub>r=1</sub></sup>&frasl;<sub>a<sub>ij</sub></sub>) - 1</sup>&frasl;<sub>v - 1</sub> [patch<sub>ij</sub>])<br>c<sub>ijr</sub> is the contiguity value for pixel r in patch ij, a<sub>ij</sub> is the area of the respective patch (number of cells) and v is the size of the filter matrix<br>Range [0,1]; contig_mn equals 0 whenever there are one-pixel patches. It reaches the max when patches are fully connected.<br>(Lagro, 1991; McGarigal et al., 2012)</td><td>Increase</td><td>2.3*10<sup>-1</sup></td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>**Fig. 5.** National forest ecosystem condition assessment through the top 10 LM-based indicators. A. Normalized feature trends over the years considering observations; B. Normalized feature trends over the years for predictions. Legend: Obs - Observed maps; Pred - Predicted maps; effective mesh size - class level (mesh_c), splitting index - class level (split_c), mean radius of gyration - landscape level (gyrate_mn_l), largest patch index - class level (lpi_c), mean of core area - class level (core_mn_c), core area percentage of landscape - class level (cpland_c), Simpson's evenness index - landscape level (siei_l), mutual information - landscape level (mutinf_l), Simpson's diversity index - landscape level (sidi_l), and mean of contiguity index - landscape level (contig_mn_l).</figcaption></figure>
fluctuates slightly from 0.5 to 0.6 overtime for the predictions. The Shannon Diversity Index, decreases slightly from 0.4 in 2000 to 0.3 in 2018, in observations and a similar downward trend, declining to 0.2 by 2054.

## 4. Discussion

The findings of this research deepen the knowledge of how forest future dynamics could be anticipated considering the spatial influence of forest/non-forest configuration, anthropogenic factors, and environmental pressures. Recognizing the most effective LM-based indicators of forest condition offers valuable insights into the modifications in forest functioning from various spatial perspectives, and impacts on biodiversity and natural capital (Skidmore et al., 2015). Additionally, developing a list of indicators that could be acquired on a cheaper and more systematic ground supports the implementation and monitoring of environmental frameworks such as the SEEA-EA and the CBD-UN. Nationally, these results provide consistent information to integrate the national ecosystem accounts, deepening the knowledge of the interactions between forest resources and human activities, and supporting decision-makers and landscape managers in predicting forest dynamics scenarios, AGB carbon stocks, and fragmentation effects.

### 4.1. Prediction of national forest distribution

Forest development scenarios allow the assessment of the capacity provision of ES, their actual and future use and values, support forest inventory, map forest habitat, and measure forest fragmentation (Zhang et al., 2022). The accuracy and quality of predictions were summarized through well-established functions and evaluated through spatial pattern analysis of forest patch centroids. Calibration and sensitivity analysis helped adjust model outputs to match observed data more closely, evaluating how changes in input data or parameter settings affected model predictions. Identifying these influential factors helped prioritize data processing efforts, refine model structure, improve prediction accuracy, and enhance model outputs' reliability. Furthermore, it is worth highlighting that the accuracy of predictions heavily relies on how well the CA model captures the essential features and dynamics of the system under study (Mengist et al., 2021), which includes the choice of cell states, neighbourhood configurations, rule selection, and boundary conditions. The choice of transition rules governing cell state changes is critical. If the rules are too simplistic or fail to capture relevant dynamics, the model's predictions may be inaccurate. Conversely, overly complex rules may lead to computational inefficiency without significantly improving predictive accuracy (Roodposhti et al., 2020). The spatial and temporal resolution of the CA grid also affects prediction accuracy. Coarse grids may overlook fine-scale patterns and interactions, while excessively fine grids may lead to computational challenges without significantly improving prediction quality. Turner et al. (1993) explored these effects on landscape patterns, revealing that as grain size increases, the landscape tends to become more uniform with high dispersion and small patch sizes. The prediction quality decreased as the difference between the start and end years increased. Furthermore, the simulations were conducted conservatively, avoiding overestimations on the predictions of forest conditions.

The fragmentation of ecological landscapes leads to more patches, smaller patch sizes, and reduced connectivity (Jaeger, 2000). Between 2000 and 2018, the number of patches increased by 22.5 %, while the mean patch area decreased by 10.3 %. However, predictions for 2036 and 2054 suggest that despite a potential decrease in the number of patches, forest connectivity may remain due to enhanced accessibility, likely from reforestation, patch consolidation, or better management. From 2000 to 2018, observed maps show rising patch numbers and mean betweenness and degree values, while mean closeness decreased, indicating growing fragmentation, but future projections suggest stable or improving connectivity. Furthermore, the connectivity of landscapes

is highly scale-dependent (Keitt et al., 1997), and the scale of the analysis can considerably impact the results, as spatial models are known to be grain and scale-sensitive (Haines-Young et al., 2012). When spatial data is aggregated, the loss of information may lead to increased homogeneity and reduced variance in coarse-resolution data, thereby modifying average values per spatial unit (Gret-Regamey et al., 2014).

Nevertheless, even anticipating an increase in forest areas influenced by initiatives, such as REDD (Reducing Emissions from Deforestation and Forest Degradation), subsidies and payment schemes, due to population growth, climate change and increased demand for natural resources including fertile lands, forest fragmentation and biomass losses are globally expected (Qasim and Csaplovics, 2023). The comparison of predictions with observed data unveiled the potential of CA systems in capturing environment-human interactions, understanding both the current trends and potential future scenarios. These assessments allow for more informed decision-making across diverse contexts, helping to maintain ecosystem functioning and resilience in an increasingly fragmented world.

### 4.2. Analysis of forest indicators and ecosystem conditions

A major finding for landscape planning and EA was that, without the inclusion of LM as a means to assess ecosystem conditions, the actual potential of modelling future landscape scenarios would be poorly realistic and expansive. This research contributed to a more systematic appraisal of the potential for predictive maps to provide long-term information to assess forest ES and the effective use of modern technologies to detect changes and their influence on natural capital and biodiversity. The natural capital and the delivery of ES are strongly affected by LULC changes, and ecological fragmentation (Jaeger, 2000). Measuring landscape characteristics and analyzing changes in the adjacency of forests and the number of patches supports understanding the consequences of these land conversions (Osborne and Alvares-Sanches, 2019). Fragmentation alters the structure of the ecological landscape and negatively affects natural capital, devaluing ES (Holzwarth et al., 2020). The fragmentation of extensive forests into smaller units with increasing area perimeter ratios indicated growing exposure to human activities (Kundu et al., 2022). Moreover, fragmented habitat reduces biodiversity and can bring a greater risk of the disappearance of some plant species and animals (Mengist et al., 2021). It can also have augmented edge effects, which can adversely influence the survival of the native species and enhance the chance of biological invasion by alien species with similar ecological niches (Adhya et al., 2022). Disturbance conditions give rise to unique landscape dynamics, highlighting their scale-dependent nature (Dou et al., 2023). Fragmentation is not only about the physical separation of landscapes but also the disruption of ecological continuity, connectivity, productivity, functional integrity, and ecological diversity (McGarigal et al., 2012).

Understanding the metrics and information that are more useful for ecosystem assessment is essential to maintaining viable ES in the long term (Meddens et al., 2022). The selected indicators provide systematic measurements of forest dynamics, revealing changes in fragmentation, connectivity, and patch structure. These indicators are critical in tracking ecosystem health, planning conservation efforts, and guiding effective forest management strategies across various contexts. Their applicability extends beyond forest ecosystems, offering insights into urbanization, agricultural expansion, and habitat restoration, where the spatial arrangement and connectivity of patches are equally critical (Arora et al., 2021). When fragmentation occurs, the edge area ratio of a forest's unit area increases, raising the probability of anthropogenic interference (Jiang et al., 2022). The largest patch index and the mean radius of gyration were the edge and area metrics that demonstrated this influence in the results. Between 2000 and 2018 the first metric decreased 52.6 %, while the second reduced by 3.6 %. The split and mesh metrics characterize the anthropogenic penetration of landscapes from a geometric point of view (Jaeger, 2000). Split increased by 28.5 %
over the years, meaning the number of forest patches also increased, while the mean mesh size decreased by about the same percentage, indicating a change in patch structure. Abundance and evenness were quantified by diversity measurements of forest/non-forest areas. Both, the Simpson's evenness index and Simpson's diversity index decreased around 28.2% over the years. Lower values are an indicator of forest extent reduction and homogenization of the landscape. The selected metrics allow for a comprehensive understanding of how forest ecosystems are being modified spatially, which affects biodiversity, habitat connectivity, and ES such as carbon storage and sequestration, soil erosion control, and water regulation and provision.

Accurate AGB estimations can significantly improve the efficiency of forest management, deepen the understanding of the forest carbon cycle, quantify and monitor terrestrial carbon stocks, and highlight landscape change trends (Han et al., 2022). The mean AGB stock at the country level was calculated with the total class area metric, which in 2000 was 2.3*10⁶ ha against 1.6*10⁶ ha in 2018, representing a loss of 39.9% over 18 years. Consequently, the national forest AGB carbon stock in 2018 was approximately 40 Mt. lower than in 2000. Carbon accounting studies value the amount of sequestered carbon over time, which helps assist policymaking regarding human-environment interactions and future adaptive strategies. They also support identifying where forests might be preserved to maintain their carbon stores and where changes in land use would minimize carbon emissions.

As forest fragmentation poses significant threats to ES and biodiversity, it is crucial to anticipate and manage its impacts. By modelling interactions between individual cells, CA models can simulate complex landscape patterns and predict spatial distribution and fragmentation over time, depicting how small-scale interactions (e.g., deforestation in one area) lead to large-scale fragmentation patterns over time. For instance, CA can predict how urban encroachment, agriculture expansion, or road development might further divide forest landscapes, exacerbating fragmentation. Such models offer spatially detailed forecasts of future landscape configurations, enabling policymakers to visualize potential fragmentation “hot spots” and focus conservation efforts where they are most needed. When combined with LM-based indicators, CA models can measure forest conditions, track ecosystem health over time, and monitor how these metrics might change in response to land-use policies or conservation actions.

Landscape development scenarios support analyzing the capacity of a landscape to provide ES and to understand the connections between natural capital, society, and decision-making processes (Karasov et al., 2020). Non-forest development is characterized by reduced natural forest areas favouring built-up, artificial, and agricultural areas. The increased competition for land use driven by the expansion of clean energy deployments, such as solar farms and wind energy installations, is also becoming a significant factor in forest fragmentation. These renewable energy projects require large, often contiguous areas of land, which can lead to conflicts with forest conservation efforts. As clean energy infrastructure grows in response to global decarbonization goals, competition for space between these deployments and natural landscapes intensifies. This dynamic presents new challenges, as the need to reduce carbon emissions through renewable energy can inadvertently contribute to land-use pressures, potentially fragmenting forests and other critical ecosystems. Balancing the land demands of clean energy with forest conservation goals requires integrated spatial planning and innovative land-use strategies. Landscape development scenarios through CA models can help simulate these scenarios, providing insights into how land competition may evolve, and guiding the identification of optimal locations that minimize environmental impacts.

Further studies should expand this approach to test LM-based indicators within forest types and consider other satellite-derived datasets with higher spatial and temporal resolution. Also, a deeper understanding of land use opportunity cost within forests, agriculture and clean energy is needed to explore the possible future development scenarios. Additionally, research efforts are required to address the impacts

on related ecosystems, including freshwater and groundwater, and assess the effects on stock and capacity to supply dependent ecosystems.

## 5. Conclusions

Sustainable development at the country scale is challenging, mainly due to land conversions based on poorly informed decision-making processes. This study evaluated the contribution of ready-to-use Copernicus data and CA models to predict future maps of forest dynamics and estimate ecosystem indicators to support national forest ecosystem physical accounting. Ten LM-based indicators of forest ecosystem condition were statistically determined through PCA and selected and ranked through feature importance assessment, plus AGB forest carbon stock (the eleventh quantified metric). The identified top 10+1 measurements deepen the knowledge of the spatial influence of forest dynamics, describing landscape heterogeneity based on its composition (amount and diversity) and configuration (spatial arrangement). Besides, changes between 2000 and 2018 were mapped for mainland Portugal, and future scenarios were predicted for 2036 and 2054. Those estimations are essential for emissions inventories, forest management, and a deepened knowledge of the carbon cycle, as well as mechanisms responsible for terrestrial sources and sinks of carbon. These advancements support more systematic estimations of carbon stock accounts by directly detecting changes in biomass.

### CRediT authorship contribution statement

**Bruna Almeida:** Writing – review & editing, Writing – original draft, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Pedro Cabral:** Writing – review & editing, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Catarina Fonseca:** Writing – review & editing, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Artur Gil:** Writing – review & editing, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Pierre Scemama:** Writing – review & editing, Methodology, Investigation, Formal analysis, Data curation, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgements

This work was supported by national funds through FCT (Fundação para a Ciência e a Tecnologia), under the project - UIDB/04152/2020 (DOI: 10.54499/UIDB/04152/2020) - Centro de Investigação em Gestão de Informação (MagIC)/NOVA IMS). We are grateful to the Institut français du Portugal for fostering Franco-Portuguese collaborations through scientific periods in France (BIC-2023), and the IUEM for hosting the first author for the main research period. We thank the anonymous reviewers for providing valuable feedback on this work.

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.scitotenv.2024.177527.

### Data availability

Data will be made available on request.
## References

Adhya, T., Bagaria, P., Dey, P., Muñoz, V.H., Weerawardana Ratnayaka, A.A., Thudugala, A., Aravind, N.A., Sanderson, J.G., 2022. On the edge: identifying priority areas for conservation of fishing cat, a threatened wetland felid, amidst rapidly altering freshwater landscapes. bioRxiv, 2021-2022.

Almeida, B., Cabral, P., 2021. Water yield modelling, sensitivity analysis and validation: a study for Portugal. ISPRS Int J Geoinf 10. https://doi.org/10.3390/ijgi10080494.

Anselmetto, N., Sibona, E.M., Meloni, F., Gagliardi, L., Bocca, M., Garbarino, M., 2022. Land use modeling predicts divergent patterns of change between upper and lower elevations in a subalpine watershed of the Alps. Ecosystems 25, 1295-1310. https://doi.org/10.1007/s10021-021-00716-7.

Arora, A., Pandey, M., Mishra, V.N., Kumar, R., Rai, P.K., Costache, R., Punia, M., Di, L., 2021. Comparative evaluation of geospatial scenario-based land change simulation models using landscape metrics. Ecol. Indic. 128. https://doi.org/10.1016/j.ecolind.2021.107810.

Bateman, I.J., Mace, G.M., 2020. The natural capital framework for sustainably efficient and equitable decision making. Nat Sustain 3, 776-783. https://doi.org/10.1038/s41893-020-0552-3.

Belo-Pereira, M., Dutra, E., Viterbo, P., 2011. Evaluation of global precipitation data sets over the Iberian Peninsula. J. Geophys. Res. Atmos. 116. https://doi.org/10.1029/2010JD015481.

Bordt, M., 2018. Discourses in ecosystem accounting: a survey of the expert community. Ecol. Econ. 144, 82-99. https://doi.org/10.1016/j.ECOLECON.2017.06.032.

Bosell, H., Krieger, H., 1991. Simulation model of natural tropical forest dynamics. Ecol. Model. 59, 37-71. https://doi.org/10.1016/0304-3800(91)90127-M.

Brandes, U., 2001. A faster algorithm for betweenness centrality*. J. Math. Sociol. 25, 163-177. https://doi.org/10.1080/0022250X.2001.9990249.

CBD-UN, 2022. Decision Adopted by the Conference of the Parties to the Convention on Biological Diversity 15/4. Kunming-Montreal Global Biodiversity Framework.

Chave, J., 1999. Study of structural, successional and spatial patterns in tropical rain forests using TROLL, a spatially explicit forest model. Ecol. Model. 124, 233-254. https://doi.org/10.1016/S0304-3800(99)00171-4.

Chrysafis, I., Korakis, G., Kyriazopoulos, A.P., Mallinis, G., 2020. Predicting Tree Species Diversity Using Geodiversity and Sentinel-2 Multi-Seasonal Spectral Information. Sustainability 2020, Vol. 12, Page 9250 12, 9250. doi: https://doi.org/10.3390/SU12219250.

Coelho, M.F., Ramos, A., de Lima, I., Trigo, R., 2013. Seasonal changes in daily precipitation extremes in mainland Portugal from 1941 to 2007. Reg. Environ. Chang. 14. https://doi.org/10.1007/s10113-013-0515-6.

Copernicus Programme. CLC 2018 — Copernicus Land Monitoring Service. https://land.copernicus.eu/pan-european/corine-land-cover/clc2018.

Csardi, G., Nepusz, T., 2006. The igraph software package for complex network research. InterJournal Complex Systems 1695.

Cunha, J., Campos, F.S., David, J., Padmanaban, R., Cabral, P., 2021. Carbon sequestration scenarios in Portugal: which way to go forward? Environ. Monit. Assess. 193, 1-14. https://doi.org/10.1007/s10661-021-09336-Z/FIGURES/7.

de Brito, H.C., Rufino, I.A.A., Djordjević, S., 2021. Cellular automata predictive model for man-made environment growth in a Brazilian semi-arid watershed. Environ. Monit. Assess. 193, 1-19. https://doi.org/10.1007/s10661-021-09108-9.

de Lima, M.I.P., Santo, F.E., Ramos, A.M., Trigo, R.M., 2015. Trends and correlations in annual extreme precipitation indices for mainland Portugal, 1941-2007. Theor. Appl. Climatol. 119, 55-75. https://doi.org/10.1007/s00704-013-1079-6.

Dou, X., Guo, H., Zhang, L., Liang, D., Zhu, Q., Liu, X., Zhou, H., Lv, Z., Liu, Y., Gou, Y., Wang, Z., 2023. Dynamic landscapes and the influence of human activities in the Yellow River Delta wetland region. Sci. Total Environ. 899, 166239. https://doi.org/10.1016/j.SITOTENV.2023.166239.

Edens, B., Maes, J., Hein, L., Obst, C., Siikamaki, J., Schenau, S., Javorsek, M., Chow, J., Chan, J.Y., Steurer, A., Alfieri, A., 2022. Establishing the SEEA ecosystem accounting as a global standard. Ecosyst. Serv. 54, 101413. https://doi.org/10.1016/j.ECOSER.2022.101413.

Elhaik, E., 2002. Principal component analyses (PCA)-based findings in population genetic studies are highly biased and must be reevaluated. Sci. Rep. 12, 14683. https://doi.org/10.1038/s41598-022-14395-4.

ESRI, 2023. ArcGIS Pro - ESRI. Environmental Systems Research Institute.

European Environment Agency, 2023. Copernicus Land Monitoring Service. https://land.copernicus.eu/en.

Farrell, C., Coleman, L., Kelly-Quinn, M., Obst, C., Eigenraam, M., Norton, D., ODonoghue, C., Kinsella, S., Delargy, O., Stout, J., 2021. Applying the system of environmental economic accounting-ecosystem accounting (SEEA-EA) framework at catchment scale to develop ecosystem extent and condition accounts. One Ecosystem 6. https://doi.org/10.3897/oneeco.e.6.65582.

Fitts, L.A., Russell, M.B., Domke, G.M., Knight, J.K., 2021. Modeling land use change and forest carbon stock changes in temperate forests in the United States. Carbon Balance Manag. 16, 20. https://doi.org/10.1186/s13021-021-00183-6.

Fleming, A., O'Grady, A.P., Stitzlein, C., Ogilvy, S., Mendham, D., Harrison, M.T., 2022. Improving acceptance of natural capital accounting in land use decision making: barriers and opportunities. Ecol. Econ. 200, 107510. https://doi.org/10.1016/j.ECOLECON.2022.107510.

Fletcher, R., Fortin, M.-J., 2018. Land-cover pattern and change. In: Fletcher, R., Fortin, M.-J. (Eds.), Spatial Ecology and Conservation Modeling: Applications with R. Springer International Publishing, Cham, pp. 55-100. https://doi.org/10.1007/978-3-030-01989-1_3.

Fonseca, F., de Figueiredo, T., Vilela, A., de Carvalho, A.L., Almeida, E., Nunes, L., 2019. Impact of tree species replacement on carbon stocks in a

Mediterranean mountain area, NE Portugal. For. Ecol. Manag. 439, 181-188. https://doi.org/10.1016/j.foreco.2019.03.002.

Foody, G.M., 2020. Explaining the unsuitability of the kappa coefficient in the assessment and comparison of the accuracy of thematic maps obtained by image classification. Remote Sens. Environ. 239, 111630. https://doi.org/10.1016/j.rse.2019.111630.

Franklin, J.F., Spies, T.A., Pelt, Van, R., Carey, A.B., Thornburgh, D.A., Berg, D.R., Lindenmayer, D.B., Harmon, M.E., Keeton, W.S., Shaw, D.C., Bible, K., Chen, J., 2002. Disturbances and structural development of natural forest ecosystems with silvicultural implications, using Douglas-fir forests as an example. For. Ecol. Manag. 155, 399-423. https://doi.org/10.1016/S0378-1127(01)00575-8.

Freeman, L.C., 1978. Centrality in social networks conceptual clarification. Soc. Networks 1, 215-239. https://doi.org/10.1016/0378-8733(78)90021-7.

Gabriel, K.R., S., R.R., 1969. A new statistical approach to geographic variation analysis. Syst. Zool. 18, 259-278.

Goetz, S.J., Baccini, A., Laporte, N.T., Johns, T., Walker, W., Kellndorfer, J., Houghton, R.A., Sun, M., 2009. Mapping and monitoring carbon stocks with satellite observations: a comparison of methods. Carbon Balance Manag. 4, 2. https://doi.org/10.1186/1750-0680-4-2.

Golub, A., Hertel, T., Lee, H.L., Rose, S., Sohngen, B., 2009. The opportunity cost of land use and the global potential for greenhouse gas mitigation in agriculture and forestry. Resour. Energy Econ. 31, 299-319. https://doi.org/10.1016/j.RESENECON.2009.04.007.

Grét-Regamey, A., Weibel, B., Bagstad, K.J., Ferrari, M., Geneletti, D., Klug, H., Schirpke, U., Tappeiner, U., 2014. On the effects of scale for ecosystem services mapping. PLoS One 9, e112601. https://doi.org/10.1371/JOURNAL.PONE.0112601.

Haines-Young, R., Potschin, M., Kienast, F., 2012. Indicators of ecosystem service potential at European scales: mapping marginal changes and trade-offs. Ecol. Indic. 21, 39-53. https://doi.org/10.1016/j.ecolind.2011.09.004.

Han, J., Hu, Z., Mao, Z., Li, G., Liu, S., Yuan, D., Guo, J., 2022. How to account for changes in carbon storage from coal mining and reclamation in eastern China? Taking Yanzhou coalfield as an example to simulate and estimate. Remote Sens. 14, 2014. https://doi.org/10.3390/rs14092014.

Harris, L.D., 1984. Island biogeography theory and the preservation of biotic diversity. University of Chicago Press, Chicago. https://doi.org/10.7208/978022619950.

Hein, L., Bagstad, K.J., Obst, C., Edens, B., Schenau, S., Castillo, G., Soulard, F., Brown, C., Driver, A., Bordt, M., Steurer, A., Harris, R., Caparrós, A., 2020. Progress in natural capital accounting for ecosystems. Science 197 (367), 514-515. https://doi.org/10.1126/science.aaz8901.

Hesselbarth, M.H.K., Sciaini, M., With, K.A., Wiegand, K., Nowosad, J., 2019. Landscapemetrics: an open-source R tool to calculate landscape metrics. Ecography. https://doi.org/10.1111/ecog.04617.

Hewitt, R.J., Compagnucci, A.B., Castellazzi, M., Dunford, R.W., Harrison, P.A., Pedde, S., Gimona, A., 2020. Impacts and Trade-Offs of Future Land Use and Land Cover Change in Scotland: Spatial Simulation Modelling of Shared Socioeconomic Pathways (SSPs) at Regional Scales.

Hewitt, R.J., Díaz-Pacheco, J., Moya-Gómez, B., 2019. A cellular automata land use model for the R software environment. https://doi.org/10.31235/OSF.IO/SUY2C.

Hewitt, R.J., Roodposhti, M.S., Bryan, B.A., 2022. There's no best model! Addressing limitations of land-use scenario modelling through multi-model ensembles. Int. J. Geogr. Inf. Sci. 36, 2352-2385. https://doi.org/10.1080/13658816.2022.2098299.

Hoffmann, J., Muro, J., Dubovyk, O., 2022. Predicting species and structural diversity of temperate forests with satellite remote sensing and deep learning. Remote Sens. 14, 1631. https://doi.org/10.3390/rs14071631.

Hojo, A., Avtar, R., Nakaji, T., Tadono, T., Takagi, K., 2023. Modeling forest above-ground biomass using freely available satellite and multisource datasets. Ecol Inform. 74. https://doi.org/10.1016/j.ecoinf.2023.101973.

Holzwarth, S., Thonfeld, F., Abdullahi, S., Asam, S., Da Ponte Canova, E., Gessner, U., Huth, J., Kraus, T., Leutner, B., Kuenzer, C., 2020. Earth observation based monitoring of forests in Germany: a review. Remote Sens. 12, 3570. https://doi.org/10.3390/rs12213570.

Houghton, R.A., Hall, F., Goetz, S.J., 2009. Importance of biomass in the global carbon cycle. Eur. J. Vasc. Endovasc. Surg. 114. https://doi.org/10.1029/2009JG000935.

ICNF, 2023. Instituto da Conservação da Natureza e das Florestas [WWW Document]. URL https://sig.icnf.pt/portal/home/ (accessed 1.24.23).

Jaeger, J.A.G., 2000. Landscape division, splitting index, and effective mesh size: new measures of landscape fragmentation. Landsc. Ecol. 15, 115-130. https://doi.org/10.1023/A:1008129329289.

Jiang, F., Sun, H., Ma, K., Fu, L., Tang, J., 2022. Improving aboveground biomass estimation of natural forests on the Tibetan plateau using spaceborne LiDAR and machine learning algorithms. Ecol. Indic. 143, 109365. https://doi.org/10.1016/j.ECOLIND.2022.109365.

Jokar Arsanjani, J., Helbich, M., Kainz, W., Darvishi Boloorani, A., 2013. Integration of logistic regression, Markov chain and cellular automata models to simulate urban expansion. Int. J. Appl. Earth Obs. Geoinf. 21, 265-275. https://doi.org/10.1016/j.jag.2011.12.014.

Karasov, O., Heremans, S., Külvik, M., Domnich, A., Chervanyov, I., 2020. On how crowdsourced data and landscape organisation metrics can facilitate the mapping of cultural ecosystem services: an Estonian case study. Land (Basel) 9, 158. https://doi.org/10.3390/LAND9050158.

Keitt, T.H., Urban, D.L., Milne, B.T., 1997. Conservation Ecology: Detecting Critical Scales in Fragmented Landscapes.

King, S., Agra, R., Zolyomi, A., Keith, H., Nicholson, E., de Lamo, X., Portela, R., Obst, C., Alam, M., Honzák, M., Valbuena, R., Nunes, P.A.L.D., Santos-Martin, F., Equihua, M., Pérez-Maqueo, O., Javorsek, M., Alfieri, A., Brown, C., 2024. Using the system of
environmental-economic accounting ecosystem accounting for policy: a case study on forest ecosystems. *Environ Sci Policy* 152, 103653. https://doi.org/10.1016/j.envsci.2023.103653.

Kundu, S., Pal, S., Mandal, I., Talukdar, S., 2022. How far damming induced wetland fragmentation and water richness change affect wetland ecosystem services? *Remote Sens Appl* 27, 100777. https://doi.org/10.1016/j.rsase.2022.100777.

Lagro, J., 1991. Assessing patch shape in landscape mosaics. *Photogramm. Eng. Remote Sens.* 57, 285-293.

Mačkiewicz, A., Rataczak, W., 1993. Principal components analysis (PCA). *Comput. Geosci.* 19, 303-342. https://doi.org/10.1016/0098-3004(93)90090-R.

Madrigal-González, J., Calatayud, J., Ballesteros-Cánovas, J.A., Escudero, A., Cayuela, L., Marqués, L., Rueda, M., Ruiz-Benito, P., Herrero, A., Aponte, C., Sagardia, R., Plumptre, A.J., Dupire, S., Espinosa, C.I., Tutubalina, O. V., Myint, M., Pataro, L., López-Sáez, J., Macía, M.J., Abegg, M., Zavala, M.A., Quesada-Román, A., Vega-Araya, M., Golubeva, E., Timokhina, Y., Bañares de Dios, G., Granzow-de la Cerda, I., Stoffel, M., 2023. Global patterns of tree density are contingent upon local determinants in the world's natural forests. *Communications biology* 2023 6:1 6, 1-6. doi:https://doi.org/10.1038/s42003-023-04419-8.

Marchette, D., 2004. Random graphs for statistical. *Pattern Recogn.* https://doi.org/10.1002/047172209X.

Marques, I.G., Nascimento, J., Cardoso, R.M., Miguéns, F., Teresa Condesso De Melo, M., Soares, P.M.M., Gouveia, C.M., Kurz Besson, C., 2019. Mapping the suitability of groundwater-dependent vegetation in a semi-arid Mediterranean area. *Hydrol. Earth Syst.* Sci. 23, 3525-3552. https://doi.org/10.5194/hess-23-3525-2019.

Matula, D., Sokal, R., 2010. Properties of Gabriel graphs relevant to geographic variation research and the clustering of points in the plane. *Geogr. Anal.* 12, 205-222. https://doi.org/10.1111/j.1538-4632.1980.tb00031.x.

Matula, D.W., Sokal, R.R., 1980. Properties of Gabriel graphs relevant to geographic variation research and the clustering of points in the plane. *Geogr. Anal.* 12, 205-222. https://doi.org/10.1111/j.1538-4632.1980.tb00031.x.

May, R.M., 1975. Species abundance patterns. *ecological diversity* 190, 1086. https://doi.org/10.1126/science.190.4219.1086.a. *Science* (1979).

McGarigal, K., Cushman, S.A., Ene, E., 2012. FRAGSTATS v4: Spatial Pattern Analysis Program for Categorical and Continuous Maps.

McGarigal, K., Marks, B.J., 1995. FRAGSTATS: spatial pattern analysis program for quantifying landscape structure. In: Pacific northwest Research Station | PNW - US For.

Meddens, A.J.H., Steen-Adams, M.M., Hudak, A.T., Mauro, F., Byassee, P.M., Strunk, J., 2022. Specifying geospatial data product characteristics for forest and fuel management applications. *Environ. Res. Lett.* 17. https://doi.org/10.1088/1748-9326/ac5ee0.

Mengist, W., Soromessa, T., Feyisa, G.L., 2021. Monitoring Afromontane forest cover loss and the associated socio-ecological drivers in Kaffa biosphere reserve, Ethiopia. *Trees, Forests and People* 6. https://doi.org/10.1016/j.tfp.2021.100161.

Moreno, N., Wang, F., Marceau, D.J., 2010. A geographic object-based approach in cellular automata modeling. *Photogramm. Eng. Remote Sens.* 76, 183-191. https://doi.org/10.14358/PERS.76.2.183.

Mpakairi, K.S., Dube, T., Dondofema, F., Dalu, T., 2022. Spatio-temporal variation of vegetation heterogeneity in groundwater dependent ecosystems within arid environments. *Ecol Inform* 69, 101667. https://doi.org/10.1016/j.ecoinf.2022.101667.

Mulatu, K.A., Mora, B., Kooistra, L., Herold, M., 2017. Biodiversity monitoring in changing tropical forests: a review of approaches and new opportunities. *Remote Sens.* https://doi.org/10.3390/rs9101059.

NASA, 2023. Earthdata | Earthdata [WWW Document]. URL https://www.earthdata.nasa.gov/ (accessed 3.2.23).

Nowosad, J., Stepinski, T.F., 2019a. Information theory as a consistent framework for quantification and classification of landscape patterns. *Landsc. Ecol.* 34, 2091-2101. https://doi.org/10.1007/s10980-019-00830-X/FIGURES/2.

Nowosad, J., Stepinski, T.F., 2019b. Information theory as a consistent framework for quantification and classification of landscape patterns. *Landsc. Ecol.* 34, 2091-2101. https://doi.org/10.1007/s10980-019-00830-x.

Osborne, P.E., Alvares-Sanches, T., 2019. Quantifying how landscape composition and configuration affect urban land surface temperatures using machine learning and neutral landscapes. *Comput. Environ. Urban. Syst.* 76, 80-90. https://doi.org/10.1016/j.compenvurbsys.2019.04.003.

Ozturk, D., 2015. Urban growth simulation of Atakum (Samsun, Turkey) using cellular automata-Markov chain and multi-layer perceptron-Markov chain models. *Remote Sens.* 7, 5918-5950. https://doi.org/10.3390/rs70505918.

Paul, K.I., Roxburgh, S.H., England, J.R., 2022. Sequestration of carbon in commercial plantations and farm forestry. *Trees, Forests and People* 9, 100284. https://doi.org/10.1016/j.TFP.2022.100284.

Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., Duchesnay, É., 2011. *Scikit-learn: machine learning in Python. J. Mach. Learn. Res.* 12, 2825-2830.

Pickard, B.R., Meentemeyer, R.K., 2019. Validating land change models based on configuration disagreement. *Comput. Environ. Urban. Syst.* 77, 101366. https://doi.org/10.1016/j.compenvurbsys.2019.101366.

OpenStreetMap Wiki [WWW Document], 2023 URL https://wiki.openstreetmap.org/wiki/Main_Page (accessed 5.22.23).

Qasim, M., Csaplovics, E., 2023. Comparative study of forest biomass and carbon stocks of Margalla Hills National Park, Pakistan. *Forest Sci Technol* 19, 139-154. https://doi.org/10.1080/21580103.2023.2208141.

R Core Team, 2021. R: A Language and Environment for Statistical Computing.

Rapinel, S., Bouzillé, J.-B., Oszwald, J., Bonis, A., 2015. Use of bi-seasonal landsat-8 imagery for mapping marshland plant community combinations at the regional scale. *Wetlands* 35 (6), 1043-1054. https://doi.org/10.1007/s13157-015-0693-8.

Rijal, S., Rimal, B., Acharya, R.P., Stork, N.E., 2021. Land use/land cover change and ecosystem services in the Bagmati River basin. Nepal. *Environ Monit Assess* 193, 651. https://doi.org/10.1007/s10661-021-09441-z.

Romme, W.H., 1982. Fire and landscape diversity in subalpine forests of Yellowstone National Park. *Ecol. Monogr.* 52, 199-221. https://doi.org/10.2307/1942611.

Roodposhi, M.S., Hewitt, R.J., Bryan, B.A., 2020. Towards automatic calibration of neighbourhood influence in cellular automata land-use models. *Comput. Environ. Urban. Syst.* 79, 101416. https://doi.org/10.1016/J.COMPVURBSYS.2019.101416.

Santos, P.L.M., Couto, F.T., Dias, S.S., Ribeiro, N.A., Salgado, R., 2023. Vegetation fuel characterization using machine learning approach over southern Portugal. *Remote Sens Appl* 101.1017. https://doi.org/10.1016/J.RSASE.2023.101017.

Savary, P., Foltête, J.-C., Moal, H., Vuidel, G., Garnier, S., 2021. graph4lg: a package for constructing and analysing graphs for landscape genetics in R. *Methods Ecol. Evol.* 12, 539-547. https://doi.org/10.1111/2041-210X.13530.

Simpson, E.H., 1949. Measurement of diversity. Nature 1949 163:4148 163, 688-688. doi:https://doi.org/10.1038/163688a0.

Skidmore, A.K., Pettorelli, N., Coops, N.C., Geller, G.N., Hansen, M., Lucas, R., Mücher, C.A., O'Connor, B., Paganini, M., Pereira, H.M., Wang, T., Wegmann, M., 2015. Environmental science: agree on biodiversity metrics to track from space. Nature 523, 403-405. https://doi.org/10.1038/523403a.

Thanh Noi, P., Kappas, M., 2017. Comparison of random Forest, k-nearest neighbor, and support vector machine classifiers for land cover classification using Sentinel-2 imagery. *Sensors* 2018, Vol. 18, page 18 18, 18. doi:https://doi.org/10.3390/S18010018.

Turner, M.G., Romme, W.H., Gardner, R.H., O'Neill, R.V., Kratz, T.K., 1993. A revised concept of landscape equilibrium: disturbance and stability on scaled landscapes. *Landsc. Ecol.* 8, 213-227. https://doi.org/10.1007/BF00125352/METRICS.

United Nations, et al., 2021. System of Environmental-Economic Accounting-Ecosystem Accounting (SEEA-EA).

Wang, Q., He, H.S., Liu, K., Zong, S., Du, H., 2023. Comparing simulated tree biomass from daily, monthly, and seasonal climate input of terrestrial ecosystem model. *Ecol. Model.* 483, 110420. https://doi.org/10.1016/J.ECOLMODEL.2023.110420.

Wilcox, B.A., Murphy, D.D., 1985. Conservation strategy: the effects of fragmentation on extinction. *Am. Nat.* 125, 879-887.

Zarandian, A., Badamfirouz, J., Musazadeh, R., Rahmati, A., Azimi, S.B., 2018. Scenario modeling for spatial-temporal change detection of carbon storage and sequestration in a forested landscape in northern Iran. *Environ. Monit. Assess.* 190. https://doi.org/10.1007/s10661-018-6845-6.

Zhang, Y., Sharma, S., Bista, M., Li, M., 2022. Characterizing changes in land cover and forest fragmentation from multitemporal Landsat observations (1993-2018) in the Dhorpatan hunting reserve. Nepal. *J For Res (Harbin)* 33, 159-170. https://doi.org/10.1007/s11676-021-01325-9.