Future Scenarios for Land Use in Chile: Identifying Drivers of Change and Impacts over Protected Area System

: Chile is a country that depends on the extraction and export of its natural resources. This phenomenon has exacerbated different processes of transformation and disturbance of natural and human ecosystems. Land use change has become a key factor for the transformation of ecosystems, causing consequences for biodiversity conservation. In this study, current and future (2030, 2050 and 2080) land use categories were evaluated. Land use projections were analysed together with models of ecosystem distribution in Chile under different climate scenarios, to ﬁnally analyse different dynamics of land use change within the protected areas system. In all the scenarios evaluated, land use projections showed an increase in the areas of industrial forest plantations and urban areas and a decrease in natural and agricultural areas could be expected. In relation to ecosystem modeling, vegetational formations located in the center and south of the country could be expected to decrease, while vegetational formations in the north and center of the country could extend their surface area. Inside Chile’s protected area network, anthropic disturbances are currently undergoing expansion, which could have consequences for ecosystems and protected areas located in the central and central–south zones of Chile.


Land Use and Land Cover Context
Humans have satisfied the need to obtain resources for their survival, such as water, wood, food, and fibre, by transforming natural ecosystems [1,2]. Approximately half of the Earth has been directly transformed by human activities [3][4][5][6]. It is estimated that, on a global scale, between 1700 and 1992, 1621 million hectares were used for agriculture, of which 885 million were native forests and 565 million ha were savanna/grassland/steppe, causing agriculture to be the most extensive land use type on the planet [1,3]. Between 1992 and 2015, there was a 115% increase in urban lands (38.0 million ha), of which 64% occupied agriculture areas, and the remaining 36% as a function of areas related to natural forests, scrublands, grasslands, and bare soil. This caused the formation of new crop areas from the conversion of natural forests (56%) and scrublands (30%), while the remaining 14% was from pastures and other uses. As a result, new cropland areas occupied a 139% larger area than those lost as a result of previous urbanisation [7].
Land use and land cover (LULC) changes have been considered an essential factor in global change [8][9][10][11][12]. Their study requires integrating social and environmental factors and disciplines [13][14][15][16][17]. Timber extraction, settlements, urbanisation, cultivation in various Land 2021, 10, 408 3 of 21 up to 44.7% of the cultivated and inhabited area of native forests and natural areas. [51]. In northern Chilean Patagonia (73 • 20 W-39 • 25 S and 71 • 59 W-41 • 14 S) in 1985, forest plantations covered 0.7% of this area and by 1998 they reached 4% of the area, reaching their maximum in 2004 (6%) [48]. The cities of Chillán and Los Ángeles (central southern Chile) have experienced a decrease in agricultural areas of −22.8% (1826 ha) and −32% (2914 ha), respectively, due to urban growth [52]. Urbanisation processes in Chile have had an impact on the configuration of ecosystems, causing problems of biodiversity loss and substitution of agricultural areas [53]. Cities in northern Chile (Antofagasta, Calama, Copiapó, and La Serena-Coquimbo conurbation) show urban growth trends during the 2005-2013 period, given their population increase related to mining factors that have attracted local and international labour [54]. In cities in central southern Chile, there is a strong association between the growth of cities and the expansion of forestry plantations, as well as the conversion of traditional crops to export plantations [52].
There is a clear relationship between anthropogenic activities, climate change, and their impact on ecosystems [2,55,56]. Climate change is an essential driver of the distribution of biomes on Earth, and several predictions indicate that in the future, climate change could alter the environmental conditions of different ecosystems in response to human activities [26]. This effect could have indirect consequences on LULC change, and at the same time, this phenomenon could cause human activities to accelerate the impact of climate change on a global scale [2,30]. LULC change could influence ecosystem processes by altering the land surface's physical properties, such as albedo, roughness, and evapotranspiration [55]. Therefore, it is essential to combine species/ecosystem correlative models with LULC change models to make projections considering different factors that may affect potential ecosystem distributions [26]. In this sense, anthropogenic activities represent one of the main factors that can lead to the extinction of species due to the disturbance and transformation of ecosystems, affecting their diversity and functioning [6].
The need to supply food and nutritional requirements and other products from agricultural and forestry systems [57,58] has resulted in more evident anthropogenic transformations of ecosystems in recent decades. For example, at least 75% of the land area that can be used for anthropogenic activities has already been disturbed [28]. Therefore, LULC change scenarios could be the most relevant tool for assessing the potential effects of LULC change on ecosystems [23,25]

Land Use Change and Ecosystem Distribution
On a global scale, major ecosystem transformations have been observed. For example, in the last decade, forest regions have lost 2.3 million square kilometres, and other regions related to natural biomes (shrublands, savannas, grassland) have lost~45% of their original extent [6]. In Latin America, since 1980, there has been an accelerated deforestation process, which has increased considerably since 1990 [59]. It is estimated that the South American tropical and subtropical native forests suffered losses of approximately 24,330 km 2 /year between 1990 and 2000 [60]. In the period from 2000-2010, global tropical and subtropical forest losses were estimated to be at least 32%, and at least half of these were in South America [24,61], which is the region with the most significant loss of native forests due to deforestation [24], a phenomenon that showed losses of woody vegetation between 20% and 80% and an increase of 359,738 km 2 in agricultural and forest plantations [58]. The expansion and intensification of agriculture have caused significant vegetation losses, concentrated mainly in tropical and subtropical forests. Given the accelerated growth of the population, it is expected that the demand for food will increase, and with that increase, there will be an expansion of agricultural areas or at least greater specialisation and intensification in the areas dedicated to food production [62].

Land Use Change and the System of Protected Areas
Chile has been subject to multiple laws and decrees that have promoted the conservation of its natural areas, and on the other hand, have encouraged forestry and natural Land 2021, 10, 408 4 of 21 resource production. In 1931, the first Forestry Law was decreed, which encouraged the national production of forest species, exempting producers from property and inheritance taxes. [46]. In 1970, the National Monuments Law was published, which has been subject to permanent modifications (the last one in 2019) and is intended to protect the cultural and natural heritage declared as National Monuments [63]. In 1974, Decree Law 701 (DL 701) was passed to promote the forestry industry through the protection of forest lands [46] and the establishment of subsidies [64] covering 75% of the establishment value of the forest plantations [46,65,66]. In 1976, the General Law of Urban Planning and Construction was published, which attempted to regulate urban growth at the expense of other land uses [67]. In 1984, the National System of State Protected Wildlife Areas (SNASPE) was created, and in 1994 Chile's National Environmental Law was published [63], encouraging the creation of protected areas on private property that are subject to the same obligations as the National System of State Protected Areas [68]. In 2008, Chile's Native Forest Law was created with the objective of creating a legal framework to protect, recover, and improve native forests [69] ( Figure S1).
In Chile, there are a total of 148 terrestrial protected areas (87 with effective management), which cover 151,000 km 2 , representing 21.3% of the country's area. The protected area network is concentrated in Chile's southernmost part, a relatively inefficient system for the country's central and northern regions [63]. Only small and scattered areas are protected [70]. The Chilean protected area network covers two of the 34 world biodiversity hotspots. These hotspots are zones recognised as having a high level of endemism and, given the high anthropogenic impact they experience, are considered global conservation priorities [71]. Therefore, on a national or regional scale, LULC change studies could generate mitigation and adaptation strategies that would empower us to effectively manage protected areas and their borders in response to possible conservation threats and opportunities.
The stability of protected areas is threatened by anthropogenic events, mainly represented by LULC changes at the borders of and within protected areas [70,72]. Protected areas and their biodiversity may also be exposed to the effects of climate change [72], thereby weakening protected areas and diminishing their ecosystem functions and stability [73].
This study aims to describe the current state of LULC in Chile, to project the future dynamics of change at a high spatial resolution (1 km) in different time steps (2030, 2050, and 2080) under four IPCC AR5 scenarios (RCP 2.6, 4.5, 6.0, and 8.5), and to relate these dynamics to the distribution of ecosystems on a national scale, characterising trajectories of future changes in both natural and anthropogenic LULC at the national scale. Finally, LULC dynamics are analysed inside the protected area network to identify potential effects on ecosystem representativeness and network integrity.

Study Area
The study area (Figure 1) comprises the continental area of Chile; the country extends over 750,000 km 2 between the coordinates 66 • E, 17 • S and 75 • E, 56 • S. For purposes of presentation and organisation of the results, the study area has been analysed in four geographic zones: North (tropical climate), Central (Mediterranean climate), Central South (temperate climate), and Austral (temperate/antiboreal climate) [74]. These zones allow us to describe LULC changes in relation to the main environmental features of the country.

LULC Modeling
An analysis of the current state of LULC change in Chile and the transformations observed in the period 1992-2011 was carried out, analysing 5 categories of land use: natural, agriculture, urban, forestry, and mining. Subsequently, LULC change was modeled for the periods 2030, 2050, and 2080, under four IPCC AR5 scenarios (RCPs 2.6, 4.5, 6.0, and 8.5) ( Figure 2). Land Change Modeler (LCM) software, available in TerrSet version 18.31, was applied; this software provides tools for evaluating and modelling LULC [23,75,76]. LCM is based on Markov chain matrices and susceptibility maps obtained by logistic regression, training machine learning algorithms, and generating prediction models [77].

LULC Data
LULC data were obtained from two different data sources, the CCI Land Cover Project [78] for 1992-2011 and the Corporación Nacional Forestal [79]; the former allowed the generation of a general dataset of LULC categories (natural, agriculture, urban, and mining), while the latter allowed the classification of LULC changes in forestry areas, which was performed for each of the administrative regions of Chile ( Figure 1) based on the years in which the forest cadastre was conducted (between 1997-2008). Finally, all the LULC categories were merged into a single raster for modelling.
The drivers of LULC change selected for the calibration process were chosen based on previous studies and represented environmental and social characteristics [32,[80][81][82][83]. Sixteen drivers were selected (Table 1), which included physiographic, climatic, and anthropic features.

Model Calibration and Validation
Three moments were identified (T 0 + T 1 = T 2 ) that allowed the model calibration to be run: T 0 corresponds to a 1992 map, which was used as the base time, and T 1 corresponds to maps of intermediate years (between 1997-2008), which vary according to when the forest cadastre was carried out by CONAF. The simulation was performed between these dates, allowing us to obtain a simulated map (T 2 ), which corresponds to the final date in the calibration process (2011); this map was subsequently compared with an observed 2011 map. The sixteen drivers selected were evaluated using Cramer's coefficient (Table 1). This coefficient is based on the chi-square statistic [23,89] and allows drivers to be selected if they contribute significantly to the explanation of the spatial distribution of the categories under study [23]. Drivers with a Cramer's coefficient greater than 0.15 better explain the change. Therefore, these drivers will have the greatest impact on the modification process and its spatial distribution [94]. We incorporated constraints in the model, and we established that the pixels over these zones (Table 1) have a value of 0, which indicates that these areas should be considered immutable, while the pixels with potential for change were assigned a value of 1, i.e., they are zones in which there is a greater probability of LULC change occurring.
Through the multilayer perceptron (MLP) method offered by LCM, transitional submodels were run (indicating the probability of changing individual LULC categories into other categories), thus generating an individual transition map that was later used for the desired prediction. The parameters used in this stage's configuration included a momentum of 0.5 and a stop order when the modelling accuracy reached 100%, the Root Mean Square(RMS) reached 0.01, or when the 10,000 interactions were completed.
LCM applied transition probability matrices that, based on the generated evaluation among the LULC maps of the years established in the input, determine the change probabilities expected for a period; therefore, these matrices were constructed by comparing the transitions generated and incorporated into the model. The model was executed for 2011, which allowed us to obtain a simulated LULC map. This simulated map was later compared with an observed map created using two methods, the receiver operating characteristic (ROC) and kappa index. ROC analysis allows the evaluation of the validity of a model that predicts the location of the presence of a class, comparing a suitability map that represents the probability of the class presence and a Boolean map (0, absence; 1, presence) of the LULC categories [95]. The kappa index evaluates the agreement between two rasters; in this case, the simulated map and the observed map were evaluated. The kappa index indicates the agreement observed between these datasets relative to the agreement expected by chance [96,97]. Kappa index coefficient values vary between 0 and 1, with the highest value indicating a good model prediction [98,99].

LULC Simulations
In the simulation of future scenarios, the variables indicated in Table 1 were considered, except for temperature and precipitation, which were changed to the variables corresponding to the four IPCC AR5 scenarios (RCPs 2.6, 4.5, 6.0, and 8.5) and for the years 2030, 2050, and 2080. These variables were used to run the transition submodels, and the probability matrices (i.e., transition probability matrices) were generated for all the years evaluated.

LULC Validation
To make an independent assessment of the LULC projections for mainland Chile, we carried out a comparative analysis with a global LULC projection. The Land-Use Harmonization (LUH2) Project [100] provides harmonised land use data for the years 1500-2100 ( Figure S2). In this dataset, we can observe historical and future projections of LULC changes at the global level at a spatial resolution of 0.5 • . Despite the differences in spatial resolution from our results, the LUH2 dataset allows us to perform visual and statistical comparisons of our models at a national scale with models at the global scale until 2080 and to make the comparison between our baseline and the LUH2 projections. To do this, averages were taken of the pixel values of the baseline period (1992-2011) and three analogous future (2080) scenarios (RCP 2.6 with SSP1, RCP 4.5 with SSP2, and RCP 7.0 with SSP3). The LULC categories included for the comparisons were those related to agricultural LULC: annual and permanent C3 and C4 annual crops as well as nitrogen-fixing crops.

LULC and the Protected Area Network
To assess LULC characterisation within protected areas, rasters of each future scenario were extracted with a mask from the protected area network using the clip tool in QGIS 3.4.4 software [101]. To quantify this change, a transition matrix was created using the "crossTabulate" function in the R library "lulcc", which allowed us to quantify and obtain information about potential LULC changes under four scenarios within the network of protected areas.

LULC and Ecosystem Modelling
Vegetational formations described in Luebert and Pliscoff (2017) [74] were used as proxies of Chilean ecosystems. Nineteen categories are described, following dominant vegetation physiognomy and macroclimate distribution ( Table 2). To model future ecosystem projections, we applied a bioclimatic modelling approach [102]. As a first step, we performed current projections of ecosystems using a climate surface dataset for the 1950-2000 period [103] and future projections (for 2030, 2050, and 2080) based on regional scenarios (RCPs 2.6, 4.5, 6.0, and 8.5) taken from CCAFS (http://ccafs-climate.org/ accessed on 3 June 2020) [104]. Future anomalies in the CCAFS dataset were matched with the current baseline following the delta method [105]. Additionally, elevation, slope, hillshade, and aspect were included as explanatory variables. All the variables have a spatial resolution of 1 km. The bioclimatic modelling step was developed in R Core Team (2018) [106] using biomod2 [107]; 19 bioclimatic variables and four topographic variables were used (Table S1). Each ecosystem was transformed to a presence map, converting the layer from a raster to a point shapefile, considering the centroid of each pixel as the point of presence/occurrence. From the total number of point occurrences, we selected a subset based on the total percentage of each ecosystem within the total study area. Pseudo-absence points (10,000 points) were randomly selected in the study area in the same proportion [108] for all ecosystems. The correlations between the bioclimatic and topographic variables were determined to eliminate multicollinearity; Pearson correlation with a threshold of 0.7 [109] and the variance inflation factor (VIF) method were implemented, where VIF values higher than 5 in the first instance or 10 in a second analysis could indicate multicollinearity between variables, which would cause prediction errors [110]. For the current ecosystems, biomod2 projections were performed, which were evaluated by partitioning the point dataset into subsets for calibration (70%) and assessment (30%). For the evaluation of the predictive power of the models, the true skill statistic (TSS > 0.8) and receiver operating characteristic (ROC) curve (ROC > 0.8) methods were used. After modelling the current ecosystem, an assembly method, including all the models, was implemented in Biomod2, which allowed us to reduce the uncertainty associated with models with low predictive power (<0.8 in the TSS and ROC). The models with the best predictive power were considered for projection to the current climate space and then projected to the future climate space for 2030, 2050, and 2080 under RCP 2.6, 4.5, 6.0, and 8.5 emission scenarios. The models used in all the ecosystem projections were random forest (RF), GBM, MAXENT, GLM, and ANN. The methodological framework and steps involved are described in Figure 2.
The current and future models of the ecosystems and LULC spatial distribution were analysed jointly. Using the crop and mask tools available in the R library "raster", we performed a spatial subtraction of the areas of nonnatural LULC for the current and future scenarios, thus obtaining information corresponding to which ecosystems are currently and could be in the future more vulnerable to possible changes in LULC due to anthropogenic actions. Figure 3 shows that the influence of the drivers (Cramer's V coefficient) varies according to the geographic zone to be analysed. The North Zone (Figure 3a) is characterised mainly by mining-related activities and shows tendencies towards urban growth, therefore, the most influential drivers in this zone are distance to mining fields (0.23) and distance to urban centers (0.20). The Central Zone (Figure 3b) is characterised by a greater influence of LULC related to forestry and agricultural production (Figure 4), consequently, the drivers with the greatest influence are distance to disturbances (0.26), distance to forestry fields (0.34), and distance to urban centers (0.27) and, in addition, temperature (0.30) and precipitation (0.29) have an important influence (Figure 3). For the Central South Zone, the drivers with the greatest influence are distance to forest fields (0.36), temperature (0.34), precipitation (0.32), distance to disturbances (0.30), and distance to the road network (0.28) (Figure 3), in this zone, there is no presence of mining activities and the most important LULCs at present and which could be more important in the future are agricultural and forest uses (Figure 4). The Austral Zone could experience few LULC transformations, characterised mainly by a slight increase in agricultural areas ( Figure 4); the most important drivers to explain land use change transformations are distance to disturbances (0.22) and distance to agricultural fields (0.23) (Figure 3).

LULC Change
Our results allow us to identify the trajectories for each LULC category. Figures 4 and 5 present the main LULC changes between 1992 and 2011, when natural and agricultural uses decreased by −0.21% and −0.20%, respectively ( Figure 5). If we express these values based on the percentage (Figure 5a) of change and ha (Figure 5b), we can observe that natural LULCs have a percentage change of −0.98%, equivalent to a decrease of 642,500 ha. These results represent the primary LULC type with the greatest loss of area. In contrast, the LULC category that lost the second most area from 1992-2011 was agricultural use, with a total percentage change of −0.20% (~607,800 ha) (Figure 5a).  As expected, the principal increases in the study area for the 1992-2011 period occurred for urban areas and forest use areas, with 83,300 ha and 1,167,000 ha, respectively. However, depending on the percentage of change, urban use showed an increase of 39.44% from its initial area, while forestry use increased by 34.76% from 1992 to 2011 (Figure 5a).

LULC Simulation
After the model calibration and validation processes, we proceeded to perform the simulation for 2030, 2050, and 2080. We established four simulation scenarios based on RCP scenarios 2.6, 4.5, 6.0, and 8.5. Figure 4 shows the changes in LULC expected for all the scenarios and four anthropic categories, which are mainly concentrated in the central and central southern areas of the country and, to a lesser extent, in the northern and southern zones.
At the national scale, we expect general decreases in natural and agricultural LULC types in the years 2030, 2050, and 2080 ( Figure 5 and Figure S3). The natural LULC will have an area between 66 million ha and 61 million ha under the RCP 2.6, 4.5, 6.0, and 8.5 scenarios; the area of this category would decrease, to cover approximately 83% of the country's surface (Figure 5c). When analysing the individual behaviour of LULC categories by geographic zone and climatic scenario, we would expect natural land use in the North Zone to decrease by an average of 1.8% in 2080 compared to 2011. In the Central and Central South Zones compared to 2011, natural land use in 2080 could be expected to decrease by 12.5% and 16.5%, respectively. Meanwhile, in the South Zone, the decrease in 2080 could represent 1.5% in relation to 2011 ( Figure 6). For agricultural LULC, the crop area will cover between 4.3 million ha and 4.6 million ha, equivalent to approximately 6% of the country's surface area. The extent of urban use areas will increase by approximately 50%, reaching between 478,900 ha and 491,800 ha, depending on the scenario ( Figure 5). The North Zone could show gradual increases in agricultural land use, which would show a greater surface area in any climate scenario in 2080 (154% more on average than in 2011). Contrary phenomena would occur in the other geographic zones, for example, in the Central Zone, a decrease compared to 2011 of 6.5% on average in 2030, 15.74% in 2050, and 30.5% in 2080 could be expected ( Figure 6). We might expect significant growth in the areas used for forestry, which would reach 1,965,100 ha if the RCP 8.5 scenario were fulfilled by 2080 (10% of the national area). If the conditions in the RCP 2.6 scenario (2080) were reached, we could expect forestry land use to increase by 5.7 million ha compared to that in 1992. For all the scenarios, it is likely that mining areas will cover 325,900 ha by the year 2080 ( Figure 5). The growth of forest LULC areas could be observed with greater specificity in the central zones where there could be increases of 3.9% on average in 2030 and 90% in the year 2080 than in 2011. In the Central South Zone of Chile in relation to 2011, for the year 2030, an average increase of 8% of the forest area could be expected, while in the year 2080, the increase could be 115% ( Figure 6). Urban land uses could show growth in all scenarios evaluated, being higher in the Central Zone, where increases could be expected between 21% (2030), 26.9% (2050), and 133% (2080). In the North Zone, urban growth could be expected to increase between 10.6% (2030) and 139% (2080). In the Central South Zone, urban area could be expected to increase by 19% on average in 2030, 38% in 2050, and 95% in 2080 ( Figure 6).
LULC transitions for the period 2011-2080 are presented in Figure S3, and we expect that between 81% and 83% of the natural LULC types will be preserved in all scenarios (Figure 5b). Furthermore, between 1.19% and 2.55% will transition towards agricultural use, between 1.88% and 2.79% could transition towards forestry uses and, potentially, 0.27% and~0.07% of natural LULC could transition towards urban and mining uses, respectively (Figure 5c).
Agriculture-related LULCs could maintain between 3.22% and 3.65% of their area, and areas related to forestry uses will maintain~3.91%. Nevertheless, agricultural land use will transition in greater proportion to forestry-related land use, which will retain between 3.72% and 3.91%. It should also be noted that forest use could change to urban (between 0.00844% and 0.01%), natural (0.06% and 0.14%), agricultural (0.18% and 0.21%), or mining (0.0004%) LULCs ( Figure S3). approximately 1515 ha (0.014%), that of agricultural LULC decreased by~2200 ha (0.02%), and that of forestry use decreased by~4500 ha (0.03%). By 2030, in the Chilean protected area network, natural LULC may increase by~14,000 ha and~20,000 ha (−0.023%) under scenarios RCP 2.6 and RCP 8.5, respectively. In 2080, the expected LULC changes show that natural coverage will decrease by approximately 191,160 ha under the RCP 2.6 scenario, and under the RCP 8.5 scenario, the decrease in natural LULC is expected to bẽ 334,770 ha by 2080. Inside Chile's protected areas, agricultural use is expected to increase 25,035 ha and 54,948 ha for scenarios RCP 2.6 and 8.5 by 2080, while by 2030 and 2050 under the RCP 8.5 scenario, agriculture use will decrease by~98,179 ha and 39,497 ha, respectively (Table S2). Urban use may increase between~3000 ha and~12,000 ha, depending on the year and scenario evaluated. Forest use would increase by~156,000 ha in the RCP 2.6 scenario and~270,127 ha in scenario 8.5, while for 2030 and 2050, these changes are expected to be~27,600 ha and~33,000 ha, respectively. Areas related to mining use, being the most stable of all the land use categories within the protected area network, could increase by~446 ha (0.003%) in all the scenarios evaluated by 2080. Table S2 shows the LULC dynamics for the 2011-2080 period for national parks, a stricter category for protection according to IUCN criteria. The main changes are projected for national parks located between the Biobío and Los Lagos regions. These changes are driving the expansion of forestry plantations.

LULC and Protected Area Network
The complete list of Chile's protected areas with LULC change dynamics is presented in the Supplementary Data (Figures S4-S18). The distribution of LULC change projections within the protected area network is shown in Figure S19. Figure S2 shows the proportion of expected change by 2080 (values range between −1 and 1, the maximum proportions of change). The expected change for RCP 4.5 is concentrated in the Central (Mediterranean climate) and Central South (temperate climate) Zones. The proportion of change was between 0.115 and 0.14 ( Figure S20). These are the two zones for which we would expect the most difference for the assessed period; these differences would be concentrated and have greater severity in scenario RCP 7.0. Of the scenarios evaluated in this research, RCP 7.0 is most similar to RCP 6.0 and RCP 8.5, with change ratio values of~0.161. The analysed LUH2 categories are mainly concentrated in the central and southern regions, where, according to our data (Figure 4), transitions to forest and agricultural uses are most expected, as well as in the northern part of the country, where the changes observed mainly indicate transitions to agricultural and mining uses. Our model shows trajectories related to a decrease in agricultural LULC areas in the study area, which contrasts with the analysed LUH2 categories, which show an increase in the area of agricultural LULC; this difference could be explained by the inclusion of different variables (socioeconomic, resource use, productivity, or energy) in each of the global scenarios of the LUH2 model.

LULC and Ecosystem Modelling
There could be latitudinal (southward) and longitudinal (towards the Andes mountain range) movements of vegetational formations in the different evaluated scenarios. These movements are more precise in the ecosystems located in the northern and central zones of Chile. In contrast, the southern ecosystems could show syndromes of redistribution within the climatic patterns most suitable for developing the species they harbour. In the northern zone, the distribution of the desert may not change excessively. Similarly, the dynamics of dunes with aerophytes may vary depending on the scenario and their distribution may shift towards higher zones, at about sea level. On the other hand, ephemeral herbaceous vegetation is one of the ecosystems that could increase in area by a greater magnitude relative to its current area ( Figure S21), with its distribution extending towards the north and south. Ecosystems such as desertic scrub and thorny forest could increase in area, occupying areas to the north and south, respectively (Figure 7 and Figure S22).
These movements are more precise in the ecosystems located in the northern and central zones of Chile. In contrast, the southern ecosystems could show syndromes of redistribution within the climatic patterns most suitable for developing the species they harbour. In the northern zone, the distribution of the desert may not change excessively. Similarly, the dynamics of dunes with aerophytes may vary depending on the scenario and their distribution may shift towards higher zones, at about sea level. On the other hand, ephemeral herbaceous vegetation is one of the ecosystems that could increase in area by a greater magnitude relative to its current area ( Figure S21), with its distribution extending towards the north and south. Ecosystems such as desertic scrub and thorny forest could increase in area, occupying areas to the north and south, respectively (Figure 7 and Figure S22).  The thorny shrubland and sclerophyllous shrubland could decrease in surface area, as a result of random behaviour by the former and movements towards the south and restrictions in the coastal zone of the latter. The surface area of the sclerophyllous forest could increase, showing variation in its extent depending on the scenario evaluated, with changes towards the south and north and in the mountainous regions ( Figures S6 and S7).
The distribution of the deciduous forest could undergo slight changes, increasing in surface area, with transitions towards the south of the country and the coastal zone. The surface area of the deciduous shrubland could increase significantly, forming new zones longitudinally throughout the country. The surface area of the broad-leaved forest could increase, with migrations towards the south of the country. The area of coniferous forest will not change without greater changes in its distribution zone (Figure 7 and Figure S21). Evergreen forest may be one of the most stable vegetation formations, as little variation is expected in its area, except under the RCP 8.5 scenario, in which the area of evergreen forest would decrease by~50% by 2080. The area of evergreen shrubland could increase, with migratory patterns towards coastal areas and fjords ( Figure S22).
LULC change in Chile could severely affect the ecosystems of the central and southern parts of the country due to the potential increases in the forestry and agricultural industries and the growth of cities. Figure 4 shows the zones that would transition to anthropogenic LULCs in the future. In this context, dunes with aerophytes would be the least impacted vegetational formation, while thorny forest and sclerophyllous forest would be the ecosystem types with the most important anthropic influence in terms of LULC change (Table S3).

Discussion
Our results show that in the period from 1992-2011 in Chile, there was an increase in forest plantations of 1.1 million ha, at a rate of approximately 75,000 ha per year. Maestripieri et al. [92], indicate that in Chile, from 1975-2007, the increase in forest plantations reached 95,000 ha/year, totalling 2.2 million ha in this period. According to Salas et al. [111], in the period 1970-1980, between 160,000 and 200,000 ha of natural forests were replaced by exotic forest plantations, representing native vegetation losses between 145,000 and 180,000 ha/year. These changes can be explained by the application of Decree 701 in the 1970s. This decree offered the application of a subsidy for new forest plantations and reforestation, with the legal guarantees that neither the lands that benefited from the subsidy nor other lands would be expropriated and the taxes levied on the existing agricultural and forestry properties would be eliminated [50,69,92,112,113]. The result of this decree was an increase in new plantations of Monterey pine (Pinus radiata) and eucalyptus (Eucalyptus globulus). According to Salas et al. [111], plantations of these species represented 62% and 31% of the total plantation area in Chile, respectively, in 2016.
This increase in forest plantations could explain the decrease in areas dedicated to agricultural uses. Agricultural use areas decreased by~607,000 ha, as indicated in the results and shown in Figure 5. These areas were replaced by forest and urban uses, a finding partially corroborated by the VII National Agricultural and Forestry Census of Chile, which could show a decreasing trend in the area of agricultural LULC in this period. This census indicated that between 1997 and 2007, there was a decrease in agricultural areas of~195,000 ha, losses that are mainly concentrated in annual and permanent crops (~100,000 ha) and permanent forage and rotation crops (~53,000 ha) and, in this period, cereal crops decreased their area by approximately 170,000 ha. The difference in the decrease in the agricultural area between the census and our model could be due to the calibration period, since the intercensal period is 1997-2007, while our model was calibrated in the period 1992-2011. In this sense, the variables used in the calibration period could also have affected, as well as the spatial resolution of the model, which could have affected the detection of smaller farms at this scale (1 km).
One of the limitations of this research was the lack of experiments performed at a similar scale in Chile. Studies have been done for specific cases, such as cities, counties, or basins, in which deforestation or the expansion of agricultural areas has been evaluated. For example, for the Biobío region, Torres et al. [51] show an intense change from agricultural land to forest plantations, with the loss of traditional crops such as wheat. In the 1997-2007 period, the agricultural area was reduced by 24.7%. Maestripieri et al. [92] show a tendency to displace traditional agricultural areas and cut native forests. These activities lead to the growth of forest areas and urban centres. Between 1999 and 2008, a rate of afforestation and reforestation equal to 9.98% per year was generated. Therefore, it is expected that by 2026, under an intensive scenario, there will be an increase in plantations of P. radiata and E. globulus, until there are~54,000 ha of plantations dedicated to timber export, replacing the native species located in this area (22,039 ha to 23,045 ha). Therefore, existing studies support the main results of this analysis.
Regarding the model projections for 2030, 2050, and 2080, given an adequate model calibration process, future change is expected, since the proportions of change are related to the change observed in the previous process. For all the RCP scenarios, an increase in forest plantations of over 50% is foreseen in its original area in place of agricultural and natural LULC. Areas dedicated to agricultural uses will decrease. However, these will shift in latitude, replacing natural LULC. In all the scenarios, urban use areas will increase by more than 55%. Henriquez Ruiz [52] indicates that in Los Ángeles and Chillán, the urban area is expected to increase 57% by 2035. For the Valparaíso region, Sandoval and Romero [114] forecast a rise in the urban area of 57% by 2025, while our results indicate an increase of 53.94% in this region; therefore, if we consider the heterogeneous growth of the urban areas in different Chilean regions by 2080, the results may seem to be in contrast with those of other studies, which addressed particular cases in greater detail and specific drivers that influence increases in specific land uses.
The effects of anthropogenic pressures around protected areas are being closely observed, e.g., projections of LULC change by 2050 and 2100 made for 50 km buffers around the region's protected areas under different RCP scenarios. Reports show that cultivated lands occupy between 7% and 8% of the spaces surrounding the protected areas in the region, and this could remain relatively stable over time [115]. In the United States, Radeloff et al. [116] show that between 2000 and 2030, approximately 17 million homes could be built in the vicinity (50 km buffer) of protected areas and wilderness areas and that this expansion could reduce natural habitats around protected areas by 12% to 16% by the year 2051 [117].
Chile's protected area network is at risk of becoming biological islands [70], reducing its functionality, effective area, and the amount of conservation space inside and outside reserves, and increasing the exposure of its borders to anthropic pressures [115]. In the context of climate change, these events could lead to more significant problems and changes in ecosystems [72,115]. Our results show that both around and within protected areas, one could expect patterns of change in LULCs, which, depending on the climate scenarios analysed, will vary on their surface. However, they will not cease to be dangerous for the stability of ecosystems.
Finally, this study, which incorporated different future scenarios, shows that a diverse set of factors is involved in LULC change processes, and these factors are difficult to interpret and incorporate into studies of this nature. The understanding of these factors, which act as driving forces, is essential to understanding LULC change dynamics. It would allow the generation of public policies that will enable us to make relevant decisions about the development of the country and the effective maintenance of Chilean protected areas.

Conclusions
This research considered different methods to achieve its objectives. For LULC change modelling, we used a model (LCM) that we considered adequate for the scale and spatial resolution. One of the limitations of this research was that there were restrictions in relation to the temporal consistency of the information, since, for the elaboration of the maps used in the calibration of the models, it was necessary to do it for each of the administrative regions of Chile for each year provided in the forest cadastre. This caused a delay in the processing of the information, however, it showed positive results in the sense that it allowed particular approximations to each region that could be analysed as an individual study area. The model calibration showed a good fit, which allowed us to meet the proposed objectives and adequately answer the research question posed. In relation to the ecosystem distribution models in Chile, 10 different models were evaluated, of which five were used, since they showed the best performance in their evaluation.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/land10040408/s1, Figure S1. Timeline of the main policies that regulate land use in Chile and promote ecosystem protection. Figure S2. The proportion of expected change by 2080 according to the Land-Use Harmonization (LUH2) Project. Figure S3. Main LULC transitions (percentage %) simulated between 1992 and 2080 for the RCP 2.6, 4.5, 6.0, and 8.5 scenarios. Figures S4-S18. LULC change (2080) within the protected area network by region in Chile. Table S1. Bioclimatic and topographic variables. Figure S19. Distribution of LULC change projections within the protected area network. Figure S20. Expected change for each category assigned according to the Land-Use Harmonization (LUH2) Project. Figure S21. Comparison of the range of changes between RCP 2.6 (by 2030) and RCP 8.5 (by 2080). Figure S22. Ecosystem area trajectory by 2030, 2050, and 2080 for RCPs 2.6, 4.5, 6.0, and 8.5. Table S2. Transitions (ha and percentage of the area) between categories by scenario in the protected area network in Chile. Table S3. Comparison of surface area (in ha) between models with and without anthropic areas due to LULC change (current scenario; 2050, RCP 6.0; and 2080, RCP 8.5).