www.mdpi.com/journal/ijgi/ Quantifying Landscape-Scale Patterns of Temperate Forests over Time by Means of Neutral Simulation Models

Several studies attempt to describe changes in the spatial patterns of forests over time, resorting to the comparison of landscape pattern indices (LPI), but new methods for quantifying landscape differences in a statistical context are necessary. In this paper, we quantified and assessed the statistical significance of the forests pattern changes, which have occurred since the end of WWII in Central Italy (Isernia). To do this; based on the proportion of forest cover (pi) and contagion (H) of three land cover maps (1954–1981–2006); we generated 100 forest maps with predictable results through the midpoint displacement algorithm. Then, for both observed and simulated maps, we computed a set of LPI (number of patches, cohesion, largest forest patch index and area weighted mean shape index) and we derived their empirical distributions; finally, we compared the empirical distributions using the non-parametric Kruskal-Wallis test. Our results show significant changes in the spatial pattern of forests and underline the process of natural forest re-growth, which, in the area, is constrained by “remnants” of traditional activities. The adopted approach could be extended to a large ensemble of landscapes and spatial scales and could become a standard procedure when comparing patterns in time.


Introduction
Landscapes change in both structure and function [1] and said dynamism is mainly driven by changes in management practices responding to social, political and economic forces [2].The analysis of temporal and spatial changes represents one of the most challenging problems in landscape ecology.Describing changes in different land cover types through time may be crucial, both for preserving biological diversity and related ecosystem services, and for developing general landscape models, which are useful in ecosystem management and environmental policies [3,4].The current state of landscape patterns in the world is mainly the result of centuries of evolution in land use [5].For instance, temperate forests in mountain and hilly landscapes of Europe are currently distributed in a mosaic, which is derived not only from centuries of extensive forest exploitation, but also from agriculture and cattle farming [6][7][8], followed by a recent, spontaneous reforestation process, which occurred after the World War II rural exodus and the consequent abandonment of traditional agricultural practices [9][10][11][12][13][14].
The quantitative, spatial and temporal analysis of natural forest re-growth in abandoned farmlands has acquired increasing relevance due to the effects of forest expansion on many important ecosystem functions [11].For example, forest expansion affects hydrological cycles and soil dynamics [15], climate [16,17] and biodiversity [11,18,19] at different scales.Several landscape metrics have been developed to quantify forest patterns in terms of space and time, most of which have been tested on grid-based thematic maps [5,[20][21][22].Studies aimed at quantifying forest distribution frequently employ landscape pattern indices (LPI) to measure changes in both forest cover and pattern.Changes over time are assessed by many authors simply through the comparison of LPI [23][24][25].Current research suggests that the comparison of raw LPI values should be avoided [26], since they are sensitive to scale [27], land cover proportions [28], spatial resolution [29,30], spatial extent [31] and land-cover misclassification [32].New methods for comparing LPI values may be useful in order to add statistical context to landscape pattern analysis [27,33].As at present it is possible to clearly identify and quantify the differences between any two given dates, the statistical significance of the observed changes remains the most important and complex challenge with which to deal [28,34].In order to enrich landscape pattern analysis with statistical significance, new methods for comparing and testing LPI values have been proposed [34][35][36].In ecological research, at the landscape scale, replication cannot be considered, and consequently, distribution, expected values and variance are not available; ecologists who want to perform statistically robust comparisons of landscapes among different maps must rely upon simulations based on computer-generated models able to reproduce an expected pattern, which shares statistical properties with an empirical pattern of interest [34,35,37].
Among the spatial models developed in landscape ecology, the neutral models (NLMs) are those able to produce an expected pattern in the absence of specific landscape processes [21,[38][39][40][41][42], providing a baseline in order to evaluate the influence of landscape heterogeneity on ecological processes and vice versa [42].The first neutral model in landscape ecology was the simple random map developed from percolation theory [39,43] and created by randomly assigning habitats to a proportion, pi, of a grid map [39].Subsequently, the evolution of NLMs followed the development of fractal methods applied to percolation maps [42].The second generation of NLMs were hierarchical, random landscapes, generated by fractal curdling [44].Both the above-mentioned NLMs assume a complete spatial independence from the habitat (or cells in a map) throughout a landscape.Since, in real-life landscapes, habitats usually show a certain degree of spatial contagion, another fractal algorithm, specifically the -midpoint displacement algorithm‖ [45] has been proposed, giving way to a third generation of NLMs, whose surfaces exhibit continuous environmental variability.Successively, many attempts to test their usefulness in various fields of research have been made.The NLMs based on the midpoint displacement algorithm have been used to assess the effect of landscape fragmentation on population distribution [46], dispersal success of communities [47], insect movement patterns [48] and biodiversity [49].They were also used to model the diffusion of contamination in aquatic environments [38], the influence of correlated spatial patterns on species coexistence of plant communities [50] and the spatial pattern of disturbance [51] at different scales [52,53].Although, initially, many NLMs were theoretically applied [54], several efforts toward modeling complex landscapes and testing their differences using neutral fractals have been made [35,55,56]; therefore, we are confident that such an approach can be applied to landscape comparison over time, while carrying out a statistical significance assessment.
Considering the aforementioned points, the present work aims to describe forest cover dynamics in the hilly landscapes of Central Italy in the last 60 years, analyzing the spatial pattern of temperate forest patches of the area surrounding a small city in a rural European setting in detail (Isernia municipality).In particular, we focus on two questions: (1) How did the spatial pattern of temperate forests change over time (1954-1981-2006)?
(2) Is the spatial pattern change of these forests significantly different over time?
To properly handle such an issue, on one hand, we used a set of landscape pattern indices to describe the forest spatial pattern dynamics and, on the other hand, we compared LPI of real and simulated landscapes to assess the statistical significance of any possible differences.In order to simulate fractal maps, we chose the midpoint displacement algorithm, because it looks very promising for modeling the natural reforestation process, which occurs in many hilly landscapes, as it is able to represent continuous autocorrelated pattern variations [40].In particular, this fractal NLM allows, through the variation of both proportion and spatial contagion, the obtaining of different levels of habitat aggregation and patchiness [57].In order to better understand the dynamics of forest fragments and to offer the basis for pinpointing specific conservation actions and forest management strategies for temperate forests in Europe, a multi-temporal analysis, which takes into account the statistical significance of differences in landscape, should be implemented.

Study Area
Isernia (Central Italy) was selected for our analysis (Figure 1), because the abandonment of traditional rural activities in this area, along with recent developments in infrastructure, appear to have caused changes in land cover, which are typical of small cities in many areas of Italy and other Mediterranean countries [58].Furthermore, the municipality of Isernia represents a pilot study area of a larger project aimed at analyzing the landscape dynamics of rural areas, which have occurred since the end of World War II in Central Italy.Isernia covers ca.68 ha and, at present, consists of a small continuous urban area (22,000 inhabitants) surrounded by a hilly landscape with agricultural, semi-natural and natural land cover types: olive groves, crop fields, complex cultivation patterns and broad-leaved forests (mainly Quercuspubescens and Quercuscerris woodlands).Although these forests are very important from a conservation point of view (-Pannonian-Balkanic turkey oak-sessile oak forests‖ according to -Habitat Directive‖ EEC 43/92), no protected areas exist in Isernia.Altitudes range from 291 m above mean sea level (a.m.s.l.) to 906 m a.m.s.l., and the climate is temperate [12].During the last half century, almost 50% of the territory has changed, mainly due to an evolution in landscape management trends: agricultural intensification in the alluvial plains and abandonment of traditional land use activities on the reliefs and less accessible areas [58].In 1955, the study area was characterized by widespread agricultural lands accompanied by natural and semi-natural vegetation patches of broadleaved forests, thickets and pastures.On the contrary, more recently, broad-leaved forest cover increased in less accessible and hilly areas, substituting shrub lands, part of pre-existing grasslands and olive groves [58].

Data Analysis
We analyzed the spatial pattern of forests and assessed the significance of the observed differences over time, following the general framework proposed by Remmel and Csillag [35], as described in Figure 2.

Data Preparation
Three land cover maps, relative to the years 1954, 1981 [58] and 2006 (1:25,000), derived from aerial photographs, were used to perform the landscape pattern analysis at the regional scale (see Acosta et al. [58]).
Due to the fact that the midpoint displacement algorithms can generate only square maps [55], we delineated a representative (approximately 20% of hilly sectors of Isernia) 256 × 256 pixel square area, which was used to extract the same geographical window from each map to be compared in time (Figure 3).Next, land cover maps were rasterized with a spatial resolution of 10 m and reclassified in two categories: broadleaved oak forests and other cover types (arable land, permanent crops, pastures, shrub and herbaceous vegetation).

NLM
Firstly, two parameters necessary for running the NLM simulations were calculated for each date (1954, 1981 and 2006): proportion of forest cover (pi) and contagion (H).pi is given by the number of forested cells in relation with the total number of cells present in the landscape and ranges between 0 (no forest cover) and 1 (whole landscape covered by forests).H, or contagion [59], describes the adjacencies of forest -cells‖ and ranges from 0, when forest distribution is maximally disaggregated (no adjacencies among cells of the same class), to 100, when landscape is totally covered by forests (only forest-to-forest adjacencies).Next, for each above mentioned map, we generated 100 maps using the public domain software Qrule [57], freely available on-line (http://www.al.umces.edu/faculty/bobgardner.html).Fractal landscapes were generated using the midpoint displacement algorithm [45], as described in With [60] and With et al. [46].In short, for each simulation, a three-dimensional fractal surface with roughness controlled by H was created by the midpoint displacement algorithm; then, every fractal surface was sectioned at the appropriate elevation to create a two-dimensional landscape map with the requisite amount of forest given by pi.Finally, in order to avoid the -salt and pepper‖ artifact commonly present in fractal maps [55] and to make simulated maps comparable with real landscapes [56], we deleted all the patches smaller than the minimum mapping unit of the original maps [61].The removal of small patches didn't significantly affect the proportion of forest cover (pi reduction < 0.5%).

LPI Calculation
To analyze the spatial pattern of forests through time, a set of landscape pattern indices on both real and simulated maps was calculated by using FRAGSTATS 4.0 [62].After recording the total number of patches of forest (NP), we focused on the other three class landscape metrics, which had been previously reported as ecologically meaningful [22,63] and have been proven to be useful in describing patch spatial structure in a forested landscape context [64][65][66][67][68][69]: largest forest patch index (LFP), area weighted mean shape index (AWMSI) and patch cohesion index (COHESION).NP and LFP were selected, because they are related to forest fragmentation [33,66,67,70], defined as the breaking up of one large forest area into many smaller patches [71].The largest forest patch index (LFP) quantifies the percentage of total landscape area comprised by the largest forest patch.As the diminution of LFP over time is one of the most effective metrics for measuring forest fragmentation [67], its increment could be an effective indicator for describing the reverse process [72,73].AWMSI [74] measures the complexity of patch shape compared to a standard shape that, in raster format, attains its minimum value (AWMSI = 1) for squares.AWMSI values increase for more irregular and elongated shapes.We chose AWMSI, because of its capacity to distinguish between big, round-shaped patches, characteristic of well-preserved forests, and small, irregular patches, which often dominate in disturbed landscapes (for a review, see Haines-Young and Chopping [64]).COHESION measures the physical connectedness of forests and is commonly used for describing habitat connectivity [75,76].In conditions of natural forest re-growth, an increment in time of landscape cohesion and connectivity could be expected [69].COHESION, which ranges between 0 and 100, is minimal when the proportion of the forested landscape decreases and becomes increasingly subdivided and less physically connected.On the other hand, COHESION increases as the proportion of the landscape covered by forests increases (see McGarigal and Marks [74] for details).

Comparing Observed and Simulated LPI
By computing the selected LPI (NP, LFP, AWSI and COHESION) on simulated landscapes, their empirical distributions (sensu Fortin et al. [77]) for each date were obtained.First, we verified the ability of the midpoint displacement algorithm to produce plausible patterns and underlying processes by generating scatterplots that include LPI derived from both real and simulated maps (Figure 4).In particular, the following LPI-based two-parameter spaces were built: NP vs. COHESION and LFP vs. AWMSI.Then, we quantified the differences of LPI among the three dates (1954, 1981 and 2006) by performing a Kruskal-Wallis test, followed by a Mann-Whitney post-hoc pairwise test.

Results and Discussion
The analysis of the 1954, 1986 and 2006 maps shows consistent changes on both the abundance and spatial distribution of forests (Figure 3).During the last 60 years, the proportion of forests increased from pi = 0.016 in 1954 to pi = 0.185 in 2006, as the spatial contagion decreased from H = 0.93 (1954) to H = 0.57 (2006).From the post-war until the end of the analyzed period, we observed the natural spread of forests, a phenomenon that characterizes the large scale dynamics of oak woodlands in many hilly and mountainous landscapes of the Italian peninsula [14,58,68,78].In Mediterranean areas, the socio-economic post-World War II changes led to a decline of the number of people involved in traditional agricultural and grazing activities [78][79][80][81][82]. Consequently, natural re-growth of forests took place in abandoned lands [58,83].
The rate of forest spread, which in the first time span (1954-1981) was of 5.9%, increased to 11% in the second one  and is in contrast with previous studies, which pointed out an overall decrease in the rate change of secondary successions over time [84][85][86], because of the occupation of the most suitable sites during the first stages of colonization.In our case, in fact, the slow and gradual cessation of farmland practices, the persistence of agricultural activities [58] and the overexploited conditions of soil in abandoned farms slowed down the establishment and growth of woody vegetation [68,86].Instead, the observed increment of forest expansion rate after 1981 is probably related to two main factors: (i) the consistent abandonment of agricultural practices occurring in Isernia during this period [58], which provided new areas to be colonized and; (ii) the presence of nuclei of forest regrowth, which ecologically facilitated landscape evolution towards a more natural condition.
In Figure 4, scatter plots showing both the position of real and simulated landscapes within the LPI-based two-parameter space (NP vs. COHESION and LFP vs. AWMSI) are reported.It is important to note that real landscapes (red symbols) fall inside the simulated landscapes scatter cloud (grey symbols).Despite the fact that many authors stress the limits of NLMs to capture the structure of real landscapes [87], the use of the midpoint displacement algorithm has provided a reliable set of simulated maps, which adequately describe the spatial pattern of forests through time (Figure 4).The fact that all the LPI calculated from real landscapes fall within the area of the respective distribution, obtained from simulated landscapes, contrasts with previous works, which state the merely theoretical value of such simulations [42,54].It is important to note that the midpoint displacement algorithm constitutes the starting point for modeling more realistic and complicated scenarios of landscape change [38,40,60].In our case, we obtained landscapes with several degrees of patchiness and spatial aggregation, by tuning the parameters that control landscape simulations (pi = proportion of forest cover and H = forest spatial autocorrelation), of the values observed in the different maps.A significant change in all LPI values was observed when analyzing the spatial pattern of forests (Kruskal-Wallis test, p < 0.01).In particular, the number of forest patches (NP) significantly increased by about 66% (3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13) from 1954 to 1981 and by about 61% (13)(14)(15)(16)(17)(18)(19)(20) from 1981 to 2006, resulting in an overall 87% (3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20) increase for the entire period .The significant rise in the number of patches, caused by the establishment of several new nuclei of young forests into abandoned farmlands [68], characterizes the former stages of natural colonization of abandoned lands in Mediterranean ecosystems [70,[78][79][80]88,89].Contemporarily, we also detected a significant increment in the percentage of landscape occupied by the largest forest patch: in 1954, the percentage of landscape occupied by the largest forest patch was 1.17% (S.E.± 0.034), while it significantly increased both in 1981 (4.47%, S.E.± 0.16) and 2006 (11.26%,S.E.± 0.39).This is mainly due to the enlargement, over time, of existing vast nuclei and their coalescence with contiguous ones [14,68,[78][79][80][88][89][90].During the advanced phases of spontaneous regeneration in Mediterranean areas, forest cover tends to evolve into a more homogeneous distribution, with a general decrease in the number of patches [14,68,90] derived from the expansion and coalescence of the numerous small pioneering forest patches into a few larger ones.In contrast with the general trend observed inside a natural reserve in Tuscany, where the protection regime allowed for a complete regeneration of forests [68], in the analyzed territory (where no conservation or protection constraints exist), only a partial forest re-growth occurred.Indeed, as evident in other Mediterranean landscapes [83][84][85][86][87][88][89][90], many -remnants‖ of the traditional activities in Molise hilly landscapes are still present and constrain the distribution and natural spread of forests.Over the three -sampling dates‖, the shape of forest patches (AWMSI) significantly increased, ranging from 1.74 (S.E.± 0.035) in 1954 to 2.48 (S.E.± 0.052) in 1981, with a further increase in 2006 (4.13, S.E.± 0.091).In terms of COHESION, there was a significant increase of the physical connection of forest patches.In particular, the COHESION index was 95.69% (S.E.± 0.14) in 1954 and increased to 97.54% (S.E.± 0.075) in 1981, reaching 97.96% (S.E.± 0.056) in 2006.As previously described in Mediterranean landscapes, the re-colonization of oak forests in Molise occurred with an increment over time of shape complexity (AWMSI) [83,91] and connectivity (COHESION) values [78,82].Immediately after the abandonment of traditional agricultural practices, the fine-grained pattern, which characterizes landscapes affected by long-established agriculture, formed by small regularly shaped and isolated patches, was progressively replaced by a more coarsely grained pattern, which results in large, irregularly shaped forest patches [83,91].

Conclusions
The observed increment in extension and the significant changes in spatial distribution of forests suggest that the analyzed area underwent an intense process of natural re-colonization, which has slowly begun after World War II and which is still in progress.The phenomenon we observed could be considered as reforestation (sensu Sitzia et al. [69]), that is, the natural reestablishment of a forested landscape on disused agricultural lands following farm abandonment [70] in regions where the potential natural vegetation (sensu Zerbe [92]) is a forest.The presence of many (NP) irregular patches (AWMSI), with increasing values of connectivity (COHESION) and patch dimension (LFP), underlines a transitional stage of forest re-growth.
Although the accurate description of the huge ecological consequences of such transformations is beyond the scope of this work, we can point out possible effects, such as an increment in true forest species [93,94], a stronger connectivity for forest vertebrates [76,95], an improvement in CO 2 sink services [96], regulation of water drainage [97] and landslide prevention [98].
The NLMs effectively modeled the pattern of forests over time at a specific spatial resolution and could be also very useful for exploring future scenarios, responding in this way to the urgent need to predict natural reforestation process [86].Simulation models could also help to better understand how natural re-growth varies in space and time and its effects on landscape function [71].In particular, predicting the pattern of forest expansion is extremely important, because of its possible effects on many key ecosystem functions [99,100].For instance, the pattern of forest expansion affects watershed services [15,101,102], biodiversity [11,18,19] and climate at different scales [16,17].
Many of the insights and conclusions obtained in this study have been facilitated by the statistical framework provided by the utilization of NLMs.In particular, the application of the midpoint displacement algorithm: (1) allows for modeling a reliable set of maps, which adequately describe the spatial pattern of forests through time and which can be directly compared with real landscapes [56]; (2) generates a set of landscape replications, which account for the most relevant information; (3) defines the landscape expectations, which allows the statistical comparison of patterns through time [29,47,48].
It is important to note that the obtained results are strongly dependent on both the specific type of landscape and the chosen spatio-temporal scale.Since, by tuning the scale of analysis, the observed patterns and the underlying processes become finer or coarser, a sound study of landscape evolution over time should include the spatio-temporal scale [103].Therefore, such an approach could be extended across a large ensemble of landscapes [37,55] and spatio-temporal scales [52,53], thus providing relevant indications as to the changes in landscape structure over time and all the ecological and cultural consequences linked to this issue [99].
Furthermore, statistical analysis could become a standard method when comparing maps [20], especially in change detection procedures, so that a sound basis for developing efficient management policies of forests could be provided.

Figure 1 .
Figure 1.Map showing the Isernia district (Italy) and the location of the study area.

Figure 2 .
Figure 2. Proposed framework for comparing and testing differences of forest pattern over time.

Figure 3 .
Figure 3. Three forest cover maps extracted from a multitemporal dataset of Isernia, Italy.Each image is 256 × 256 pixels, with a spatial resolution of 10 m.The binary classification separates oak forest (dark green) from other cover types (arable land, permanent crops, pastures, shrub and herbaceous vegetation-soft grey).pi indicates the proportion of forest cover; H is the autocorrelation of the map.(a) = 1954; (b) = 1981; (c) = 2006.

Figure 4 .
Figure 4. Scatter plots among two paired combinations of observed landscape pattern indices (LPI) values for simulated landscapes (N = 100; 1954, 1981 and 2006): (a) total number of patches of forest (NP) vs. COHESION and (b) LFP vs. weighted mean shape index (AWMSI).Asterisks indicate the mean values, and red symbols show the position of real landscapes within the empirical distribution derived from simulated landscapes (in grey scale).