Simulation model of land cover changes in a post-socialist peripheral rural area: Požega-Slavonia County, Croatia

Simulacijsko modeliranje korisna je metoda za objašnjavanje i projekciju promjena načina korištenja zemljišta i zemljišnog pokrova, koje su u osnovi povezane sa širim društvenim i prirodnim elementima. Značajne promjene zemljišnog pokrova, pretežito napuštanje zemljišta, dogodile su se u perifernim ruralnim područjima diljem Središnje i Istočne Europe nakon raspada socijalizma. Jedno je takvo područje Požeško-slavonska županija u Hrvatskoj, za koju su modelirane promjene zemljišnog pokrova od 1985. do 2027. godine. Prevladavajuća promjena od 1985. do 2013. bila je sekundarna sukcesija. Simulacijski model budućih promjena zemljišnog pokrova temeljen je na ćelijskim automatima i umjetnim neuronskim mrežama i realiziran u dodatku MOLUSCE za QGIS. Validacija testnih modela pokazala je da su oni točnije predviđali kvantitetu (bolje od lokacija) promjena te da je naknadna kombinacija individualnih rezultata simulacija nešto poboljšala njihovo podudaranje s referentnim podatcima. Konačni je model predvidio da bi šume mogle pokrivati više od 60 % cijelog područja do 2027. ako se trenutni negativni demografski i ekonomski trendovi nastave u budućnosti.

Ključne riječi: umjetne neuronske mreže, ćelijski automati, promjene načina korištenja zemljišta i zemljišnog pokrova, periferno područje, simulacijski model, MOLUSCE Simulation modelling is a useful method for the explanation and prediction of land use and land cover changes, which are fundamentally interconnected with broader social and natural elements. Significant land cover changes, predominantly land abandonment, occurred in peripheral rural areas across Central and Eastern Europe after the collapse of socialism. One such area is Požega-Slavonia County in Croatia, for which we modelled land cover changes from 1985 to 2027. The prevailing change from 1985 to 2013 was secondary succession. The simulation model of future land cover changes was based on cellular automata and artificial neural networks, and was implemented in the MOLUSCE plugin for QGIS. Validation of the test models showed that they more successfully predicted the quantity (rather than the location) of the changes, and the post-simulation combination of individual results slightly improved their agreement with reference data. The final model predicted that forests would cover over 60% of the whole area by 2027 if the current negative demographic and economic trends continue into the future.

Introduction
Land use and land cover changes (LULC) result from the complex interaction of many natural and social factors (Pijanowski et al., 2002). Understanding the causes and consequences of LULC has become critical because they impact biodiversity, carbon dynamics, climate, hydrology, and livelihoods (Memarian et al., 2012).
Global trends of LULC can be divided into two opposite categories: intensification (agricultural expansion, urbanisation, deforestation, etc.) and extensification (land abandonment, afforestation, etc.) (Díaz et al., 2011). Polarisation of these trends is the most common on agricultural land, which represents the largest terrestrial biome. There is a consensus that agricultural intensification generates land degradation and reduces the quality and quantity of the services that ecosystems provide to humanity. However, the abandonment or extensification of agriculture yields two-sided consequences (Rey Benayas et al., 2007). Negative consequences of extensification include the disappearance of traditional farming practices, long-term loss of habitats of high ecological value, higher probability of wildfires, invasion of exotic species, food insecurity (Díaz et al., 2011), reduction of biodiversity, loss of landscape identity, etc. (Ruskule et al., 2012). Some of the positive consequences are improved hydrological regulation, soil recovery and erosion mitigation, higher water quality, increased fertility, fungal biomass and decomposer activity, and carbon sequestration (Rey Benayas et al., 2007).
These consequences are not always relevant in all parts of the world or are only relevant at small scales (Rey Benayas et al., 2007). Nowadays, the biodiversity of semi-natural agricultural land in Europe is almost as valuable as wild biodiversity. A central issue of Europe's traditional cultural landscapes is their instability, i.e. their dependence on a medium degree of human impact. If land use is abandoned, traditional landscapes will be displaced by spontaneous secondary succession. Conversely, human impact that is too intensive will lead to the conversion of traditional landscapes into more simplified landscapes (Plieninger et al., 2006). Since the mid-twentieth century, European peripheral predvidjeli nastavak toga trenda u idućim desetljećima (Navarro i Pereira, 2012;Ruskule i dr., 2012).
The interest in researching European peripheral areas is also visible in numerous rural typologies. Newer examples include typologies from Croatia (Lukić, 2012), the Czech Republic (Klufová, 2016), Hungary (Beluszky and Sikos, 2008), Moldavia (Tudora, 2009), Poland (Bański and Mazur, 2016), Russia (Mikhaylova et al., 2015), and Slovenia (Cosier et al., 2014). They indicate that rural areas with different developmental obstacles occupy a significant part of national territory. Although peripheral rural areas are certainly not a unique feature of Central, Eastern, and Southeast Europe, their specificity is their composition: not only remote or areas with natural disadvantages but also the so-called internal rural peripheries (TA, 2011). This unique type of countryside is described as an area with poor accessibility and scarcity of more developed urban central settlements. Many of them are characterised by weak local economies and lack of job opportunities for young people. Negative demographic processes progress, particularly out-migration and ageing of the population (TA, 2011). As a consequence, land is usually abandoned, resulting in secondary succession and afforestation.
Glavni ciljevi ovog istraživanja bili su: (1) utvrditi povijesne promjene načina korištenja zemljišta i zemljišnog pokrova i stvoriti njihov stohastični simulacijski model za jednu postsocijalističku perifernu ruralnu regiju -studiju slučaja Požeško-slavonske županije od 1985. do 2027., (2) utvrditi povezanost između kompleksne ruralne tipologije i promjena načina korištenja zemljišta i zemljišnog pokrova u području studije slučaja, (3) utvrditi pogodnost umjetnih neuronskih mreža i ćelijskih automata za modeliranje kompleksnih i nelinearnih promjena načina korištenja zemljišta i zemljišnog pokrova i (4) ispitati potencijal njihova simulacijskog modela za izradu razvojnih scenarija u području istraživanja. LULC in post-socialist Central and Eastern European countries has been an interesting research topic during the last few years. After the collapse of the socialism in these countries, one of the main LULC trends was decrease in agricultural land use due to extensification, land abandonment and afforestation from one side, or suburbanisation on the other (Václavík and Rogan, 2009;Baumann et al., 2012;Alcantara et al., 2013;Griffiths et al., 2013). Similar studies have been conducted in the area of northern Croatia (Horvat, 2013;Cvitanović et al., 2016;Jogun et al., 2017), where similar trends and driving forces of LULC were detected, like in the previously mentioned papers. The methodology of those papers was based on analysis of historical LULC, while this study aims to extend them using predictive simulation models of LULC. Namely, simulation models are one of the most effective tools to study LULC, and they can additionally support land use planning and policy (Verburg et al., 2004). Many simulation models have been created during the two last decades, under initiatives like the Land Use and Land Cover Change Project of the International Geosphere-Biosphere Program (IGBP), and International Human Dimensions Program (IHDP) (Han et al., 2015). Most of the empirical studies employing LULC simulation models have dealt with dynamic areas that have undergone significant LULC such as deforestation, expansion of built-up areas and pressure on agricultural land due to urbanisation, population growth, and economic development in the last few decades (Pontius et al., 2008). This research will offer new insights on the performance of LULC simulation model in a postwar and post-socialist rural area, which can be applied to similar studies in the future.
The main goals of the study were to: (1) detect historical LULC and create a stochastic simulation model of LULC for one post-socialist peripheral rural region: a case study of Požega-Slavonia County from 1985 to 2027; (2) ascertain the relationship between complex rural typology and LULC in the case-study area; (3) assess suitability of artificial neural networks and cellular automata for modelling complex and nonlinear LULC; and (4) examine the potential of the LULC simulation model for the creation of development scenarios in the research area.

Research area and data
Požega-Slavonia County, as a marginal rural area with diverse landscape features (lowlands vs mountains), and turbulent social conditions in the recent decades (fall of socialism, war, transition, and EU entrance) is a fine example for modelling land cover changes in similar areas of Central and Eastern Europe. It is located in the western part of the historical region of Slavonia, in eastern Croatia (Fig. 1). With a surface of 1823 km 2 (CBS, 2015) and composite development index of 32.8% of the national average (Pejnović and Kordej-De Villa, 2015), it is one of the smallest and the most underdeveloped Croatian counties.
Požega-Slavonia County had 78,034 inhabitants in 2011, a respective general population density of 43 persons per square kilometre, which was less in comparison to 76 persons per square kilometre in Croatia as a whole. The number of settlements was large (277), as a consequence of traditional agrarian exploitation, and historical and natural conditions. Most of the settlements had populations less than 200, and 15 of them were unpopulated. The share of the urban population was 41.6%, which was also below the Croatian average of 56%, confirming the rural characteristics of the area (CBS, 2011). The Požega Basin has a very long history of human inhabitation, and until the second half of the 20 th century, it had been attractive for many immigrants. However, since the 1971 population census, this area has been losing population. The western part of the county was particularly struck by the Croatian War of Independence (1991)(1992)(1993)(1994)(1995), the negative demographic and economic repercussions of which are still being felt today. Moreover, the negative trends of ageing and depopulation are set to undermine the demographic situation in the future (Njegač, 2012).
The temporal frame of the research was influenced by data availability, research methodology, and deep changes in the political and economic system that occurred in the research area. The base year for the test model was 1985, the launch year was 1999, and the target year for simulation was 2013. The final model had its base in 1999, the launch year was 2013, and the target year was 2027.
The main data source was multispectral satellite images acquired by the Landsat mission, which were downloaded from the USGS EarthExplorer service. The criteria for their selection were: coverage of the same research area; availability of all bands (Level 1 Data Product); cloud-free scene; the same season (summer) for phenological consistency; and equal and sufficient periods between images. The earliest images were acquired by Landsat 5 TM on August

Klasifikacija zemljišnog pokrova
Klasifikacijska shema sastojala se od pet klasa: voda; izgrađeno; usjevi i tlo; trava i šikara; šuma. Scene su klasificirane nadziranom metodom u dodatku Semi-Automatic Classification Plugin za QGIS, uporabom algoritma maximum likelihood. Uzorci za treniranje izabrani su na temelju tzv. true colour (crveni, zeleni i plavi kanal) i false colour RGB kompozita (kanali iz infracrvenog dijela spektra) te iz točnijih dodatnih podataka. S obzirom na to da se klasa izgrađenog zemljišta obično miješa s drugim klasama jer obuhvaća mozaik 10 th , 1985; the second set was acquired by Landsat 7 ETM+ on August 9 th , 1999; and the most recent images were acquired by Landsat 8 OLI TIRS on August 7 th , 2013. The images were downloaded in GeoTIFF format in 30 m spatial resolution, already georeferenced and georectified in the WGS 84 UTM 33N coordinate system. Thermal and panchromatic bands were not used. The digital terrain model, from which natural variables for the simulation model were derived, was retrieved from the European Environmental Agency. Its original spatial resolution was 25 m and its projection ETRS89 LAEA, so it was resampled to 30 m by the nearest neighbour method and re-projected to WGS 84 UTM 33N. The data on areas suspected to be mined were retrieved from the Croatian Mine Action Centre. Administrative boundaries were obtained from the Digital Atlas of the Republic of Croatia.

Land cover classification
The classification scheme consisted of five classes: water; built-up; crops and soil; grass and shrubs; and forest. The scenes were classified using the supervised method in Semi-Automatic Classification Plugin for QGIS, using the maximum likelihood algorithm. Training samples for classes were chosen on the basis of "true colour" (red, green, and blue bands) and "false colour" RGB composites (bands from the infrared and visible parts of the spectrum), and with more accurate ancillary data. Regarding the fact that the class of built-up land usually gets the most mixed up with HRVATSKI GEOGRAFSKI GLASNIK 81/1, 31−59 (2019.) zgrada, vegetacije i tla, klasifikacija je provedena odvojeno za izgrađeno i ostalo zemljište, prema preporuci u Horvat (2013). U tom je postupku ručno izdvojeno izgrađeno zemljište iz dodatnih podataka. Nakon konačnog klasificiranja načinjeno je modalno filtriranje s prozorom 3 × 3 ćelije kako bi se smanjio šum.
ANN je oblik strojnog učenja koji se sve više rabi zbog napretka računala i povećane dostupnosti fleksibilnih softvera (Pijanowski i dr., 2002). ANN je jedan od najuspješnijih sustava u daljinskim istraživanjima koji se primjenjuje za klasifikaciju zemljišnog pokrova (Gašparović i Jogun, 2018) i modeliranje promjena načina korištenja zemljišta i zemljišnog pokrova. Vrlo se dobro nosi s netočnim i nepotpunim podatcima te kompleksnim nelinearnim obilježjima the other classes because it incorporates a mosaic of buildings, vegetation, and soil, the classification was conducted separately for built-up and for other land, following recommendations from Horvat (2013). In this procedure, the built-up land was manually delineated from ancillary data. After finishing the entire classification, modal filtering with a 3 × 3-pixel neighbourhood size was carried out to reduce noise.
Reference samples for the classification accuracy assessment were determined independently from the classification results and training samples, by the same method used for the training samples. The quantitative measures of accuracy were omission, commission, and overall agreement, instead of a deficient Kappa coefficient (Pontius and Millones, 2011).

Simulation of land cover changes
The simulation model of LULC in this research was conducted in MOLUSCE (Modules for Land Use Change Evaluation), which is an open-source plugin for QGIS (MOLUSCE, 2018). The plugin implements the following functionalities: takes input data of land cover categories and explanatory variables; trains LULC model using input data and well-known algorithms; predicts future LULC based on the training input; and performs validation of results based on reference data from the past (GIS-Lab, 2018).
For modelling LULC transition potential, four algorithms can be used, i.e. artificial neural network (ANN), logistic regression, weight of evidence, and multi-criteria evaluation (MOLUSCE, 2018). We selected results from ANN multi-layer perceptron because they had the highest accuracy compared to the other algorithms.
ANNs are a form of machine learning that is used increasingly because of advances in computing performance and the increased availability of flexible software (Pijanowski et al., 2002). It is one of the most successful systems in remote sensing used for land cover classification (Gašparović and Jogun, 2018), and modelling LULC. It is very good at coping with incorrect and poor data, and capturing complex non-linear features in modelling processes (Li and Yeh, 2002 (Li i Yeh, 2002). ANN se sastoji od slojeva i neurona koji oponašaju strukturu ljudskog mozga i njegovu sposobnost sortiranja uzoraka i učenja iz pokušaja i pogrešaka uočavajući tako veze u podatcima (Li i Yeh, 2002). Najrašireniji oblik ANN-a u praksi je multi-layer perceptron (MLP), koji su opisali Rumelhart i dr. (1986).
Korak simulacije izvodi se ćelijskim automatima (CA), koji su najpoznatiji koncept LULC-a. Osnovna je ideja ćelijskih automata da se promjene načina korištenja zemljišta i zemljišnog pokrova mogu objasni-of layers and neurons that imitate the structure of the human brain and its ability to sort patterns and learn from trials and errors, thus observing relationships in data (Pijanowski et al., 2002). The most widely used form of ANN is the multi-layer perceptron (MLP), described by Rumelhart et al. (1986).
The MOLUSCE plugin uses classical realisation of an MLP ANN with backpropagation. It performs initial preprocessing of the input data (dummy coding of categorical variables into a set of independent variables, and normalisation of factor variables), sampling, and training (GIS-Lab, 2018). For each model, there were 1000 randomly selected samples based on 1-size Moore neighbourhood (3 × 3 window size).
In the backpropagation algorithm, momentum is used for the learning procedure, and weights corrections are calculated as: where w is a vector of neuron weights, dw is a vector of weights changes, n is an iteration number, r is learning rate, m is momentum. Training set is divided into learning set (80% of samples) and validation set (20% of samples). Learning is stochastic, where a random sample is selected from the learning set, and the weights of the network are updated during forward/backward propagation. A fitting error for a sample is calculated as: where E is a sample error, t i is the target value of a output neuron for given sample, o i is the real output value of the neuron, d is the count of output neurons (GIS-Lab, 2018). Models in this research were trained with learning rate and momentum parameters set to 0.001 in order to stabilise the learning graph. A number of iterations was set to 200 in order to avoid problems with model over-fitting.
Simulation step is performed with Cellular automata (CA), which is the most well-known LULC concept. The main idea of CA is that LULC can be explained by the current state of a "cell" and its HRVATSKI GEOGRAFSKI GLASNIK 81/1, 31−59 (2019.) ti trenutnim stanjem "ćelije" i promjenama u njezinu susjedstvu (van Schrojenstein Lantman i dr., 2011).
Kalibracija i validacija bitni su elementi simulacijskog modela. Jasna distinkcija između ovih pojmova načinjena je u Pontius i dr. (2004). Najjednostavnija metoda validacije je tzv. nul-model. Simulacijski model kalibrira se s podatcima iz vremena t 1 , a validacija se temelji na usporedbi simulirane karte za t 2 s opaženom kartom u t 2 . Nul-model pretpostavlja perzistenciju, tj. da nije došlo do promjena između t 1 i t 2 . Ako je nul-model (karta iz t 1 ) sličniji opaženoj karti iz t 2 nego simulirana karta, to znači da simulacija nije zadovoljavajuća. Ovdje je provedena trostruka usporedba karata koja podrazumijeva slaganje između parova triju karata: referentne karte za t 1 , referentne karte za t 2 i simulirane karte za t 2 . Ta trostruka usporedba karata omogućila je izdvajanje ćelija koje su bile točne zbog perzistencije nasuprot ćelijama koje su bile točne zbog promjena. Matematički izrazi za mjere definirane kombinacijom informacija o kvantiteti i lokaciji mogu se pronaći u prethodnim istraživanjima, tj. ukupno kvantitativno neslaganje i slaganje (Pontius i Millones, 2011), ukupni pomak i razmjena (Pontius i Santacruz, 2014), i faktor kakvoće projekcije (Pontius i dr., 2008). The CA implemented in MOLUSCE takes as input the next data: initial state raster (current land cover categories); factor rasters (explanatory variables); and model (transition potentials that are the output of an algorithm, e.g. an ANN). The simulator takes transition probabilities from the transition matrix, calls a model, and feeds it with initial state and factor rasters. The model calculates transition potentials of every transition class, and the simulator constructs a raster of the most probable transitions. For every transition class in the raster of the most probable transitions, the simulator searches a needed count of pixels with the greatest certainty and changes the category of the pixels (GIS-Lab, 2018). Simulations in this research were performed in two iterations, meaning that the described procedure was repeated with the first simulated map as the initial state map for the second iteration. Because the model integrates CA and ANN, it will be referenced as CA-ANN in the rest of the paper.
Calibration and validation are essential elements of the simulation model. A clear distinction between these terms was made in Pontius et al. (2004). The simplest method for validation is a null model. The simulation model is calibrated with data from time t 1 , and validation is based on the comparison of the simulated map for t 2 with the observed map in t 2 . The null model assumes persistence, i.e. no change between t 1 and t 2 . If the null model (map from t 1 ) is more similar to the observed map from t 2 than the simulated map, it means that the simulation is not satisfactory. We conducted a three-way map comparison, which considered the agreement between two pairs of three maps: the reference map of t 1 ; the reference map of t 2 ; and the simulated map for t 2 . This three-map comparison allowed us to distinguish the pixels that were correct due to persistence, versus the pixels that were correct due to change. Mathematical expressions for measurements defined by a combination of information of quantity and location can be found in previous papers, i.e. overall quantity disagreement and agreement (Pontius and Millones, 2011), overall shift and exchange (Pontius and Santacruz, 2014), and figure of merit for prediction (Pontius et al., 2008).
Considering that the most problematic class (built-up) was delineated manually, the overall classification accuracy was high: 98.7%; 95%; and 96. 9%-for 1985, 1999, and 2013 respectively. The largest errors of omission and commission were in the class of grass and shrubs (Tab. 1).
The largest change from 1985 to 1999 was a decrease in crops and soil, and an increase in grass and shrubs, and forest. From 1999 to 2013, the most signif- Najkompaktnija područja zarastanja obrađenih površina u travu i šikaru bila su u gorskim dijelovima Psunja i Pakračke gore te pakračko-lipičkog kraja, koji su bili najteže pogođeni Do-icant change was a reduction in grass and shrubs, and an increase in forest cover. The changes in built-up and water, the smallest classes, were not significant (Tab. 2).
The most compact areas of change from cultivated land to grass and shrubs were in the mountainous parts of Psunj and Pakračka gora, and the Pakrac-Lipik area, which were the hardest struck by Legend: 0 -no change, 1 -water, 2 -built-up, 3 -crops and soil, 4 -grass and shrubs, 5 -forest Legenda: 0 -bez promjene, 1 -voda, 2 -izgrađeno, 3 -usjevi i tlo, 4 -trava i šikara, 5 -šuma both the Croatian War of Independence and subsequent depopulation. Most of the extensive changes were grouped in the elevation zone of 200-400 m, while the highest mountainous parts (covered with forest, unpopulated) and the lowest parts of lowlands (agricultural areas, demographically relatively stable) were less prone to change ( Fig. 3; Fig. 4).

Land cover change variables
The initial variables for LULC in the test simulation model 1985-1999-2013 were: elevation; slope; aspect; distance from rivers; distance from lakes and reservoirs; distance from precedent LULC; general population density by settlements  0-14/population 65+) by settlements; distance from built-up settlements; distance from settlements with population greater than 1000; distance from traffic network; and distance from areas suspected to be mined.
Because of multicollinearity or low explanatory power, the following were excluded from the initial variables: aspect; population growth index; distance from traffic network; and distance from built-up settlements.
Different combinations of parameters and variables produced dozens of simulated land cover maps for 2013. Because of the stochastic nature of the model, three maps produced with optimal parameters were combined by a modal function into one final land cover map for 2013 (Fig. 5).

47
increased their overall agreement to 83.5%. The simulation models had slightly better results than the null model, but their overall validation results were similar. Much bigger differences between the models' performance were found in the prediction of quantity and location of land cover classes. Namely, all simulation models better predicted quantities of classes than the null model, while the null model better predicted locations of classes than the simulation models (Fig. 6).
It was more important to assess the accuracy of simulated changes than overall land cover. A comparison of net changes between observed and simulated maps revealed that the simulation model correctly predicted the tendency of changes in all classes except crops and soil. It predicted the largest net change especially well-an increase of forest cover (observed = 6.3%, simulated = 5.3%).

49
The validation map shows the spatial distribution of the agreement and disagreement components of the simulated LULC (Fig. 8). The most accurate predictions of LULC were in southern parts of Pakračka gora, where they were the largest and relatively linear due to the total land abandonment after the war.

Final simulation of land cover changes (2013-2027)
Transition probabilities for the period of 2013-2027 were computed with the Markov chain analysis from the period of 1999-2013 (Tab. 5).
Model calibration was done with the parameters ascertained as optimal in the test model, while dynamic social variables were changed (e.g. population density from 2011). After running the CA-ANN simulation, three of the produced maps for 2027 were selected and combined using the modal function, as in the test model (Fig. 9).
Regarding the fact that validation showed that the simulation models better predicted quantity rather than the location of LULC, the quantitative analysis likely deserves more attention (Tab. 6).

51
tion in less favourable areas and intensification in Požega Basin, where the natural (flat terrain, fertile soil) and social conditions (relatively stable population, better transport connectivity) are better for agriculture.
Allocation of simulated LULC in the period of 2013-2027 was similar to that of the previous period (Fig. 10). Therefore, the location-validity of the model could be considered satisfactory at the coarse spatial level.

Evaluation of methods and techniques
Hybrid classification of images acquired by remote sensing was used for detection of land cover. The high overall accuracy achieved (over 95%) confirmed that this is one of the most accurate classification methods (Horning, 2004).
Conversely, validation of simulation models depends highly on the accuracy of reference data, and comparisons between different models are relative (Pontius et al., 2008). Furthermore, many studies are lacking validation entirely (Pontius et al., 2004;Verburg et al., 2004;Pontius et al., 2008). Although the comparison to a null model is the simplest validation method, it is very rigid, because it is sometimes an insuperable barrier (Pontius et al., 2004). Pontius et al. (2008) compared 13 simulation model applications, and 6 of them performed worse than the respective null model. The figure of merit in this study was 20.2%, which is the average in comparison to the other results in Pontius et al. (2008), and better than that shown in more recent studies (Memarian et al., 2012;Sloan and Pelletier, 2012). However, Pontius et al. (2008) found that figure of merit of the LULC model depends on the amount of LULC, where models observing greater LULC have a greater figure of merit. This study confirmed that claim because LULC were observed in 17% of the area, which is also average.
Regarding quantity and location disagreement, results showed that the simulation models better predicted quantity rather than the location of LULC, which corresponds to previous studies (Memarian et al., 2012).
Besides the amount of LULC, their trends also impact model performance. If the trends are linear (e.g. deforestation in developing countries), the model will perform better. Here, non-linear changes due to time lags in the process of secondary succession, which are hard to incorporate into a model, were problematic. This research area is unique in comparison to previous studies, in the manner that there was an observed decrease
The number of land cover classes significantly impacts the accuracy of the simulation models. With the introduction of every new class, the probability of error exponentially increases. Hence, dichotomous models (urban/non-urban, forest/non-forest) perform better than ones using multiple classes (Li and Yeh, 2002). In this study, there were five land cover classes, which were enough to determine the most relevant processes of extensification (crops-grass-forest) and intensification in the study area. This number of classes was average in comparison to previous studies (Pontius et al., 2008).
Scientists are constantly trying to find ways to improve the results of the validation of LULC simulation models. A popular method is a multi-resolution validation, in which a model is validated at different spatial resolutions. Predicted maps at coarser resolution have a higher agreement with reference maps than at finer resolution (Pontius et al., 2004). However, this might indicate that the model is scale-dependent, rather than more accurate. We made a combination of simulation results to improve the overall accuracy because similar approaches in classification have yielded promising results (Chen, 2008). Considering that the model is stochastic, combination increases the probability of a correct simulation (e.g. if two models predicted the persistence of class A, and a third model predicted change from A to B, the combined result would retain persistence, which is more probable). The improvement due to the combination was not statistically significant, and the figure of merit of the combined map was even smaller than in some specific maps. This can be explained by components of agreement and disagreement because specific maps better predicted changes, but they did a worse job of predicting persistence than the combined map. Combination improved the overall visual balance of the maps. In the future research, it would be interesting to include pan-sharpening methods during preprocessing of input images and examine its effects on classification and LULC model accuracy (Gašparović and Jogun, 2018). istraživanja jer je opaženo smanjenje izgrađenog zemljišta, što se smatra vrlo neobičnom promjenom (Verburg i Overmars, 2007).

Evaluation of results
Today, post-socialist Croatia has similar problems with land abandonment in peripheral rural areas to those of other countries from the former Eastern Bloc, e.g. Latvia (Ruskule et al., 2012), and Slovakia (Némethová et al., 2014).
In the mountainous parts of Požega-Slavonia County, which were used for low-intensity agriculture 30 years ago, the dominant process was secondary succession, which will probably result in total afforestation by 2030. Similar processes of land abandonment were documented in previous studies in Croatia (Cvitanović, 2014), in which they were explained by social restructuring, i.e. deagrarianisation and industrialisation. Afforestation in Požega-Slavonia County in the period of 1985-2013 was influenced by depopulation, unfavourable natural conditions, transport isolation, deep political and economic changes (collapse of socialism, transition, EU accession), war, unsuccessful privatisation, etc.; which corresponds to more recent studies (Ruskule et al., 2012;Griffiths et al., 2013;Cvitanović et al., 2016).
The relation between afforestation/share of forest cover and types of settlements demonstrated that LULC could be used as an indicator of finescale development trends. Depopulation is one of the most important driving forces of land use extensification and overall development problems in Požega-Slavonia County. Spatial plans were recently updated with acts regarding land as an unrenewable resource, but new data on LULC trends and population dynamics are lacking. Accurate input data in development strategies and documents are the first step in their successful realisation, so they should be improved. This paper contains a lot of data on LULC that can be used in new spatial plans and strategies of Požega-Slavonia County.
Active strategies of land use and land cover management could be intensification, extensification, and afforestation. Sometimes, however, the most beneficial option is "re-wilding", i.e. passive management of ecological succession with the aim to restore natural ecosystem processes and reduce human control of landscapes (Navarro and Pereira, 2012). In any case, the decision on land use and
Aktivne strategije upravljanja načinom korištenja zemljišta i zemljišnim pokrovom mogu biti intenzifikacija, ekstenzifikacija i aforestacija. Međutim, katkad je najisplativije pasivno prepuštanje prirodnom tijeku razvoja, tzv. renaturalizaciji (engl. re-wilding) kako bi se obnovili prirodni procesi u ekosustavu i smanjio ljudski utjecaj na pejzaž (Navarro i Pereira, 2012). U svakom slučaju, odluka o strategiji razvoja korištenja HRVATSKI GEOGRAFSKI GLASNIK 81/1, 31−59 (2019.) land cover strategy should be based on regional studies (Rey Benayas et al., 2007). Remote sensing and GIS have had a great role in such studies, but they should not be seen as a panacea (Farrow and Winograd, 2001), nor are geographers and GIS-analysts the only competent experts in of LULC studies. An interdisciplinary, holistic approach is necessary, whereby simple visualisation of results should be the common language between all stakeholders (Herrmann and Osinski, 1999). The results of this study are meant to be used by the community, to help decide on the most appropriate land use management strategy.

Conclusion
This paper presented the application of integrated CA-ANN simulation model for analysis and prediction of land cover changes (MOLUSCE) in a post-socialist peripheral rural area.
The hybrid method of supervised classification of Landsat satellite images appeared to be very reliable, with a minimum overall accuracy of 95%. Through comparison of classified maps, significant LULC in Požega-Slavonia County from 1985 to 2013 were observed. The class of forest had the highest increase in overall landscape, from 47% in 1985 to 55% in 2013. The class of soil and crops had the greatest decrease, from 35% in 1985 to 26% in 2013. Thus, the main trend was secondary succession, characterised by the change from agricultural land to grass and shrubs in the period of 1985-1999, and a change from grass and shrubs to forest in the period of 1999-2013.

57
The integrated CA-ANN model was proven to be one of the most efficient methods for simulating complex and non-linear LULC, which are specific for the research area. The simulated maps had a better agreement than the null model in comparison to observed maps, especially in predicting the quantity of LULC. The post-simulation combination of individual results slightly improved the overall agreement of simulated land cover with the observed, but not significantly. The main weakness of the model is that it is focused on pattern resemblance and acts like a "black box", so it would be appropriate to employ some kind of agent-based model in a similar area and compare their performance.
The simulation model predicted that forest will cover 63% of Požega-Slavonia County by 2027 and that the mountainous zone between the Požega Basin and the Pakrac-Lipik area will be almost completely afforested. The potential dangers and opportunities of such LULC have been presented, therefore enabling the creation of management plans and scenarios. For example, in areas of Požega-Slavonia County that would become covered with forest, construction zones can be restricted to protect nature and let the land to become "wild" again; or they can be expanded, to stimulate immigration, along with other planning measures (e.g. subsidies for ecological agriculture).
Deep development problems in Požega-Slavonia County require focused and multi-sectoral strategy in which the realisation of spatial and regional development plans has an important role. Optimistically, this study would be useful in envisioning evidence-based policy and it is offered for evaluation, to all interested groups, in the hope that it will contribute to the welfare of society and nature in Požega-Slavonia County. Also, the study identified problems of LULC simulation modelling and management, which could be applied in other peripheral rural areas with similar development backgrounds.
This work was supported by the Croatian Science Foundation under Grant 4513 (Croatian Rural Areas: scenario-based approach to discuss planning and development, CRORURIS).