A Step-by-Step Guide to Initialize and Calibrate Landscape Models: A Case Study in the Mediterranean Mountains

Suárez-Muñoz, María; Mina, Marco; Salazar, Pablo C.; Navarro-Cerrillo, Rafael M.; Quero, José L.; Bonet-García, Francisco J.

doi:10.3389/fevo.2021.653393

ORIGINAL RESEARCH article

Front. Ecol. Evol., 23 April 2021

Sec. Biogeography and Macroecology

Volume 9 - 2021 | https://doi.org/10.3389/fevo.2021.653393

This article is part of the Research Topic Using Landscape Simulation Models to Help Balance Conflicting Goals in Changing Forests View all 14 articles

A Step-by-Step Guide to Initialize and Calibrate Landscape Models: A Case Study in the Mediterranean Mountains

Updated

A correction has been applied to this article in:

Corrigendum: A Step-By-Step Guide to Initialize and Calibrate Landscape Models: A Case Study in the Mediterranean Mountains
1. Read correction

$\r\nMaría Surez-Muoz*$ María Suárez-Muñoz^1*

Marco Mina²

Pablo C. Salazar³

Rafael M. Navarro-Cerrillo⁴

José L. Quero⁵

Francisco J. Bonet-García³

¹Department of Ecology, University of Granada, Granada, Spain
²Centre for Forest Research (CEF), Université du Québec à Montréal (UQAM), Montréal, QC, Canada
³Department of Ecology, Faculty of Sciences, University of Córdoba, Córdoba, Spain
⁴Silviculture and Climate Change Lab, Faculty of Forest Engineering, ERSAF Research Group, University of Córdoba, Córdoba, Spain
⁵Forest Systems Ecophysiology Lab, Department of Forest Engineering, Andalusian Institute of Earth System Research, University of Córdoba, Córdoba, Spain

The use of spatially interactive forest landscape models has increased in recent years. These models are valuable tools to assess our knowledge about the functioning and provisioning of ecosystems as well as essential allies when predicting future changes. However, developing the necessary inputs and preparing them for research studies require substantial initial investments in terms of time. Although model initialization and calibration often take the largest amount of modelers’ efforts, such processes are rarely reported thoroughly in application studies. Our study documents the process of calibrating and setting up an ecophysiologically based forest landscape model (LANDIS-II with PnET-Succession) in a biogeographical region where such a model has never been applied to date (southwestern Mediterranean mountains in Europe). We describe the methodological process necessary to produce the required spatial inputs expressing initial vegetation and site conditions. We test model behaviour on single-cell simulations and calibrate species parameters using local biomass estimations and literature information. Finally, we test how different initialization data—with and without shrub communities—influence the simulation of forest dynamics by applying the calibrated model at landscape level. Combination of plot-level data with vegetation maps allowed us to generate a detailed map of initial tree and shrub communities. Single-cell simulations revealed that the model was able to reproduce realistic biomass estimates and competitive effects for different forest types included in the landscape, as well as plausible monthly growth patterns of species growing in Mediterranean mountains. Our results highlight the importance of considering shrub communities in forest landscape models, as they influence the temporal dynamics of tree species. Besides, our results show that, in the absence of natural disturbances, harvesting or climate change, landscape-level simulations projected a general increase of biomass of several species over the next decades but with distinct spatio-temporal patterns due to competitive effects and landscape heterogeneity. Providing a step-by-step workflow to initialize and calibrate a forest landscape model, our study encourages new users to use such tools in forestry and climate change applications. Thus, we advocate for documenting initialization processes in a transparent and reproducible manner in forest landscape modelling.

Introduction

Forests are indispensable ecosystems for human societies. Due to their role as climate regulators, soil protectors and biodiversity hotspots, forests provide a multitude of ecosystem services and are fundamental elements in the world’s economy (Krieger, 2001; Martínez Pastur et al., 2018). The potential adverse impacts of global change on forest ecosystems emphasizes the need to understand how to manage them in the future (Lindner et al., 2014; Hof et al., 2017; Mina et al., 2017).

In recent years, the use of computational models has been increasing in forest ecosystem research (Gustafson et al., 2017; Shifley et al., 2017). Although empirical studies are of fundamental importance for process understanding, simulation models are nowadays recognized as useful tools to assess our knowledge about the functioning of ecosystems as well as essential allies when predicting future changes (Seidl, 2017). Over the past decades, a large range of models were developed to describe future dynamics in forest ecosystems (Keane et al., 2015), from stand-scale empirical simulators to more complex process-based models operating at landscape scale (Fontes et al., 2010). Because of computational constraints, models integrating fine-resolution processes (e.g., photosynthesis, specific growth functions) at large scales in a spatially explicit framework were rare, and smaller grain processes were often strongly simplified (Elkin et al., 2012). However, these constraints are constantly being reduced by the increase in computational power, allowing for the flourishing of Forest Landscape Models (FLMs) which integrate physiologically based processes from stand to landscape level (Seidl et al., 2012; De Bruijn et al., 2014; Shifley et al., 2017; Petter et al., 2020).

According to Jorgensen and Fath (2011), ecological models comprise five elements: state variables, external variables, parameters, mathematical equations, and universal constants. The mathematical equations and universal constants are implicit within the model structure, while the initial conditions of the state variables (e.g., species biomass, species age), external variables and parameter values are usually provided as inputs by the user for each simulation study. In the case of FLMs, they represent forests across the landscape in a spatially explicit way. The landscape is depicted as a set of cells for which a series of state and external variables are defined. These variables are used to define the ecological processes taking place at cell level (e.g., growth, mortality among others) and at landscape level (e.g., seed dispersal, fire spread). For comprehensive reviews on the development, structure and recent applications of FLMs see Shifley et al. (2017); Keane et al. (2015), and He et al. (2017).

The above-mentioned structural elements are essential to set up a simulation in a specific landscape. This requires the user to obtain, prepare and organise comprehensive datasets to address the two first key steps in applying FLMs: model initialization (initial conditions of state and external variables) and calibration of model parameters. The initial conditions of the state variables describe the ecosystem at the beginning of the simulation. In turn, external variables are those forces affecting the ecosystem without being internal parts of it (Jorgensen and Fath, 2011). Most FLMs require initial values of at least certain state and external variables to start a simulation. As an example, the FLM LandClim requires elevation or browsing intensity for initialization (Petter et al., 2020). In addition to biophysical conditions (e.g., soil types, climate maps or regions), essential initial conditions for FLMs are vegetation maps describing which species are present in the landscape at the beginning of the simulation time. Decisions regarding the inclusion or exclusion of certain species can be highly relevant in certain ecosystems (e.g., shrubs in the Mediterranean area). Thus, these vegetation maps are a keystone within these experiments since forest dynamics and properties (biomass, available light, regeneration, etc.) are highly driven by initial conditions and structure (Duveneck et al., 2015; Scheller and Swanson, 2015). For example, most FLMs require information of tree species and their age classes across the landscape. This information can be very challenging to obtain and estimate for large spatial scales without necessarily combining multiple and complex datasets (Zald et al., 2014). The generation of input data for FLMs can therefore require significant time and skills, and often demands complementarity with experimental research from long-term field studies (Shifley et al., 2017; Scheller, 2018).

The calibration of model parameters has been defined as one of the greatest challenges in modelling under environmental changes (Keane et al., 2015; Scheller, 2018). Model parameters are values used in model equations which represent processes (Jorgensen and Fath, 2011). Most models simulating the succession dynamics of vegetation require parameters describing the behaviour of the species present in the landscape. These parameters may differ for each model, but commonly refer to species growth characteristics, fruit and seed dispersal, reproduction strategies and absolute or relative measures of tolerance to stress factors (Huber et al., 2018). A broad range of sources can be used to fulfil these parameters, ranging from empirical case-specific data collected by the modeller to values of standard variables stored in global databases. In either case, an evaluation of model outputs to identify appropriate parameter values is usually required to ensure that the model produces plausible outcomes at the local scale (Gutiérrez et al., 2016; Duveneck et al., 2017). This evaluation of parameter values is known as calibration (Mulligan and Wainwright, 2013). During this phase, the different model sensitivity to some parameters over others should be considered (McKenzie et al., 2019).

Successional processes and long-term projections in FLM are highly sensitive to initial conditions and model parameters (Scheller, 2018). Estimation of initial condition and calibration procedures are typically described in the method section in literature studies (e.g., Scheller et al., 2005; Boulanger et al., 2017), but often not on a level of detail to allow full reproducibility or with enough information to help non-modellers to setup a new landscape from scratch. Even fewer provide access to inputs, outputs and scripts in public repositories. The aims of this study are twofold. First, documenting the process needed to initialize and calibrate a FLM step-by-step, as an example for analogous uses. Fulfilling this goal would encourage the application of FLMs as scientific tools to assess future forest dynamics and management adaptation under global change. Second, we aim at testing different initialization data – with and without shrub communities – by assessing the model ability to project forest landscape dynamics in a biome where the model has not been applied so far (Euro-Mediterranean region).

Materials and Methods

Model Description

In this manuscript, we chose LANDIS-II as our reference FLM. LANDIS-II¹ (Scheller et al., 2007) is a FLM designed to simulate forest dynamics at multiple spatial and temporal scales. It allows a wide degree of complexity depending on a set of extensions which can optionally be activated to simulate different processes such as succession, disturbances (fire, wind, herbivory, and pests) and management at different degrees of complexity (e.g., areas- and species-specific harvesting regimes, post-harvesting planting). The spatial scale (i.e., cell resolution) is defined by the user, which makes it very flexible and adaptable to a wide variety of simulation experiments. In LANDIS-II, the landscape is divided into ecoregions, which are subregions sharing similar climatic conditions and soil characteristics. Trees in each cell are represented as species-age cohorts, increasing the computational efficiency of the model (De Bruijn et al., 2014).

Particularly, we used the PnET-Succession extension v.3.4 (De Bruijn et al., 2014; Gustafson et al., 2015). This extension embeds the PnET-II ecophysiological model equations (Aber and Federer, 1992). PnET-Succession simulates forest succession in a more mechanistic way than previous approaches, representing an advantage for experiments where novel conditions such as climate change are being explored (Gustafson, 2013). In PnET-Succession, age is used to calculate cohort’s biomass at the onset of the simulation (i.e., model spin-up) and cohorts with higher biomass are given priority access to light and water (Gustafson and Miranda, 2019). Cohort biomass is assumed to be homogeneously distributed in the cell and therefore shade conditions are also homogeneous within a cell (Scheller et al., 2007). Potential net photosynthesis rate is calculated as a linear function of foliar nitrogen (FolN) and biomass growth is a result of environmental conditions such as temperature, precipitation, photosynthetically active radiation (PAR), CO₂ concentrations and—optionally—ozone concentrations (De Bruijn et al., 2014). Biomass allocation depends on compartments turnover and fraction parameters. Mortality can occur at any time if carbon reserves become limiting (non-structural carbon <1%) or when age approaches species longevity (De Bruijn et al., 2014).

PnET-Succession requires a series of generic, ecoregion- and species-specific parameters. Although many default values have been made available by the model developers and in past application studies, most parameters require calibration according to the biogeographical location of the target landscape and the tree species included in the simulations (McKenzie et al., 2019; Mina et al., 2021).

Study Area

The simulated area considered in this study is located in the south eastern part of Iberian Peninsula and it covers approximately 390,000 ha (37.2° N, 3.1° W, Figure 1). The topography is mostly mountainous, including three mountain ranges. In the southern part of the study area, Sierra Nevada spreads from east to west and contains the highest peak in the Iberian Peninsula, Mulhacén (3,478 m). In the northern part of the study area, Sierra de Arana is located in the west, while Sierra de Baza-Filabres is in the east. More than half of the study area is under protection, either as National or Natural Park, and therefore a variety of exploitation and management regimes can be found in the study area.

FIGURE 1

Figure 1. Study area location (A), orthophoto of the study area (B) and pictures of representative forest types: pine plantations (C), mixed open forest (D), oak-dominated stand (E). Shaded area in panel (B) delimitates pine plantations.

Several bioclimatic zones are found within the study area (Rivas Martínez, 1983; REDIAM, 2018). The supramediterranean zone (mean annual temperature 8–13°C) is the one covering most of the area, at low altitudes of Sierra Nevada and connecting with Sierra de Arana. An important proportion of the Sierra de Baza surface is also represented by this bioclimatic zone. Supramediterranean areas are mostly covered by a mosaic of conifer, mixed forests and sclerophyll shrubs. The highest areas represent the oromediterranean zone (mean annual temperature 4–8°C) and are covered by conifers, shrubs and sparse vegetation, except for the very high altitudes in Sierra Nevada, which encompass the cryoromediterranean zone (i.e., alpine tundra, mean annual temperature <4°C). These peaks are covered by scarce vegetation adapted to extreme climatic conditions. The rest of the study area, at altitudes commonly below 1,000 m, is mostly covered by mesomediterranean (mean annual temperature 13–17°C) and thermomediterranean zones (mean annual temperature 17–19°C). The precipitation follows a strong seasonal pattern, with dry summers and precipitation concentrated in a small number of events. Rainfall is the most common form of precipitation. Besides, snowfall at high altitudes is very important since slow melt down and subsequent infiltration into soil increases water availability for plants throughout the spring and summer season. Aspect also determines water availability due to different precipitation evapotranspiration patterns.

The study area is covered by diverse natural vegetation patches in combination with agriculture and Pinus forest plantations. Pine plantations are the dominant land use type, covering around 20% of the study area, with a minor presence of natural pine forests. These plantations were mainly established between the 1950s and 1970s as means to halt soil erosion in recently abandoned agricultural areas. The main species are Maritime pine (Pinus pinaster Aiton), Aleppo pine (Pinus halepensis Miller.), black pine (Pinus nigra Arnold.) and Scotch pine (Pinus sylvestris L.) (Bonet García et al., 2009; Pemán García et al., 2017; Mesa Garrido, 2019). Pines were planted in high densities to drastically reduce soil loss. Afterwards, favourable climate conditions and lack of appropriate post-planting management have resulted in highly dense monospecific even-age stands. As a result, these forest plantations are nowadays under extreme risk by climate change and forest pests, which has resulted in decline and massive mortality processes (Sánchez-Salguero et al., 2010, 2012a).

Almost 40% of the study area is covered with shrublands and abandoned crops with sparse natural vegetation. Some of these areas host sparse trees (mainly Quercus ilex L.), which can be highly relevant seed sources at a local scale. Moreover, in a context of climate change and further abandonment of mountain agriculture activities, these sparsely vegetated systems can be highly important to understand the succession dynamics in pine plantations for two reasons: (1) Due to climate change, currently forested areas could suffer a decline and be replaced by shrublands as these areas become less suitable to sustain high levels of biomass; and (2) Tree species could expand to shrubland areas currently dedicated to marginal activities (mountain agriculture, livestock, fuelwood, and charcoal exploitation, etc.).

Model Initialization

In this section we describe the workflow followed to produce the necessary inputs required by LANDIS-II with PnET-Succession. We first focus on the generation of the initial vegetation conditions, followed by the methodological process to build biophysical inputs. The different sources of information used in this process are listed in Table 1.

TABLE 1

Table 1. Information sources used in this study.

Initial Vegetation Conditions

Most FLMs require estimates of initial conditions in the form of vegetation maps reporting the presence of tree species at the beginning of the simulations. Concretely, LANDIS-II requires a spatial representation of species cohorts by age classes (also called “initial communities”). To generate such vegetation maps, different information sources are often combined following a complex workflow that, if not exhaustively described, is often unreproducible. Even though such workflows can be model- and site-specific, three methodological steps can be defined: (1) Plot-level information (e.g., from national or regional inventories or permanent growth plots) is processed to extract tree measurements such as basal area, age or height; (2) Polygon-level information is processed to select stand-scale, spatially explicit variables which can be linked to plot-level information such as forest type, mean age, canopy cover, etc.; (3) A spatial imputation method is applied to produce a continuous map by assigning plot-level information to polygon-level information. In the following sections, we describe the methodological details for each of the three methodological steps of the workflow, which is summarized in Figure 2.

FIGURE 2

Figure 2. Schematic view describing the workflow of the initial vegetation map generation. Tree species presence and age were extracted from plots of the Third Spanish National Forest Inventory (NFI) during the plot-level information processing phase. During the polygon-level information processing, the Spanish Forestry Map (FM) was used to classify polygons as active/inactive and no forest/forest. The Andalusian Vegetation Map (VM) was used to extract the presence of shrub communities and sparse trees. At the spatial imputation phase, plot- and polygon-level information was combined to generate a continuous map of vegetation.

Plot-level information processing

Plot-level information is necessary to select the most common tree species and some of their demographic features in the area of interest. In our case, the Third Spanish National Forest Inventory (NFI) was used for this purpose (Table 1). The NFI contains homogeneous information about forest covers in Spain by reporting data collected in a systematic network of permanent sampling plots (Alberdi et al., 2017). The plots are evenly distributed on a 1 km² grid throughout the territory and contain plot- and tree-level information for each survey period. In addition to single-tree data (e.g., species, diameter, height, form, and health status), the plot is described in terms of the three most dominant tree species contributing to canopy cover. We selected the tree species to include in the model simulations based on the total coverage value of the species within the study area.

Tree age is often not available at single-tree level and its estimation is challenging as several factors influence growth rate of individual trees, leading to very different tree characteristics for the same age. Nevertheless, LANDIS-II uses cohort age-classes as a proxy for biomass, and therefore an estimated age is required for each species across the landscape. NFI provides estimated stand age for plots within even-aged stands. In our area, even-aged stands are composed by Pinus spp. Since these plots also have associated individual tree measurements (e.g., diameter, height), we calculated the average diameter per species and plot and joined it to the assigned stand age. In order to have as many observations as possible, we expanded the considered dataset to all plots within even-age stands from surrounding regions (provinces of Granada and Almeria). Since no estimated age was available for plots within uneven-age stands (mainly Quercus spp. and Populus nigra), a semi-quantitative method was applied. We used yield tables available from the literature (Teobaldelli et al., 2010). Then, we generated an age assignment table containing correspondence rules between tree age, diameter and height for each species under consideration. This table was validated based on expertise and observations in the study area. Finally, the age assignment table was used to attribute an age class to the species in each NFI plot (see Supplementary Table 1 for more details).

Polygon-level information processing

Many government and private forestry organizations utilize cartographic products to support forest management. Polygon-level information usually contains forest variables at a landscape scale such as forest type, coverage or stand development stage. Here, we used the Spanish Forestry Map (FM) and the Andalusian Vegetation Map (VM) as the main source of spatially explicit landscape information (Table 1). The FM is a vectorial file generated by photointerpretation. It contains a series of attributes describing the forest vegetation for each polygon: polygon identifiers and surface; province and region; forest and land use characteristics (vegetation type, structural type, land use categories); name, coverage and state of the three most dominant species in dominance order; and tree and total coverage fraction. Based on these attributes, polygons were classified as active (those containing or potentially containing vegetation units useful for our purpose) or inactive (those where natural succession is hampered by human activities: crops, infrastructures, firebreaks, etc.). Active polygons were in turn classified as forest and non-forest. Forest polygons have an average size of 51 ha and a maximum of 357 ha.

The VM is also a vectorial file that contains an extensive list of attributes (140) describing the presence and characteristics of the tree species, shrubs and pastures present in each polygon (e.g., vegetation community, canopy coverage) at a higher resolution than the FM. Due to the importance of sparse trees and shrublands in our study area (see section “Study Area”), we used the VM to consider the occurrence of sparse trees in non-forest polygons in the initial vegetation map. Moreover, since the VM also provides information about shrubs, grasslands and pastures, we used a variable termed “life form” attribute to analyse the presence of shrubs within the study area. The life form attribute is based on the classification proposed by Raunkiaer (1934) and subsequently revised by Ellenberg and Mueller-Dombois (1967). According to this classification, plants can be: microphanerophytes (“evergreen perennial plants that grow between 2 and 5 m, or whose shoots do not die back periodically to that height limit”), nanophanerophytes (“evergreen perennial plants that grow below 2 m tall, or whose shoots do not die back periodically to that height limit”), chamaephytes (“evergreen perennial plants whose mature branch or shoot system remains perennially within 25–50 cm above ground surface, or plants that grow taller than 25–50 cm, but whose shoots die back periodically to that height limit”), hemicryptophytes (“perennial plants with periodic shoot reduction to a remnant shoot system that lies relatively flat on the ground surface”), geophytes (“perennial plants with periodic reduction of the complete shoot system to storage organs that are imbedded in the soil”), therophytes (“annual plants whose shoot and root system dies after seed production and which complete their whole life cycle within 1 year”). Shrub communities’ size was extracted from this classification. Hemicryptophytes, geophytes and therophytes life forms were not considered for the analysis as they mostly refer to species with short life cycles. As a result, the presence of tall (2–5 m), medium (0.5–2 m), short (<0.5 m) shrub communities was extracted for each polygon in the study area. Each of these shrub communities was parameterized individually in PnET-Succession (see details below).

Spatial imputation

Spatial imputation is applied to combine plot-level to polygon-level information. To generate the initial vegetation map suitable for LANDIS-II, we combined plot- and tree-level information (NFI) with polygon-level information (FM and VM). The final aim of this step is to produce a forest composition map containing the species-age assemblage (species and age of each cohort) in every cell within the study area (see Duveneck et al., 2015 for the description of a similar approach in North America).

First, we performed a spatial join between the forest polygons from the FM and the plots from the NFI (Figure 2). FM polygons containing one single NFI plot were assigned the species-age assemblage of the surveyed plot. FM polygons with more than one NFI plot were assigned the species-age assemblage which results from merging all plots species-age assemblages. This only occurred in a small proportion of cases: 80 out of 789 forest polygons had two NFI plots in them, 3 polygons had 3 NFI plots in them and one polygon had 4 plots in it.

Forestry Map polygons without NFI plots were analysed to identify an equivalent polygon among those which have one or more intersecting NFI plots. Polygon equivalency was analysed at three levels: (i) Full equivalency: polygons matching vegetation type, species dominance, species development state, and total coverage fraction; (ii) Partial equivalency: polygons matching vegetation type, species dominance, and total coverage fraction; and (iii) Species equivalency: polygons matching species dominance and total coverage fraction. FM polygons were assigned the species-age assemblage of their equivalent polygon. If there were more than one equivalent polygon, the polygon was assigned the species-age assemblage resulting from the merge of all possible ones. FM polygons with no equivalent polygons were further analysed based on their species composition, without considering species dominance order. These polygons were assigned the species-age assemblage corresponding to the most common species-age assemblage containing the same species as the considered polygon in the whole study area. FM polygons containing species assemblages not occurring in the previously analysed polygons were assigned one cohort of each of the species present in it. The age of this cohort corresponded to the most common age for each species among all previously analysed polygons. This procedure resulted in the description of the forest trees within the study area.

Second, polygons labelled as non-forest by FM, were imputed to include sparse trees (Figure 2). We used VM to gather information for those polygons. Since VM does not report any variable describing the age of the species, VM polygons were imputed a species-age assemblage containing the reported species with the most common age for that species in the rest of the study area. By doing so, the description of the sparse trees within the study area was completed.

Third, we imputed shrub communities within the landscape by assigning the corresponding shrub community-age assemblage to each VM polygon (Figure 2). Shrubs were not allocated to different age classes but instead we arbitrarily assigned age 10 for all polygons. Since in LANDIS-II shrub biomass should not exceed tree biomass (see section “Model Description”), the role of shrubs in the model is mainly to compete for light and water in the understory (e.g., affecting establishment), thus we believe age class assignation for shrubs was not necessary.

Finally, we obtained the initial vegetation map by combining forest trees, sparse trees and shrubs communities. The map was rasterized at 100 m resolution (1 ha cells). Each cell was labelled with a code associated to a list of unique species-cohort assemblages.

Biophysical Inputs

The FLM LANDIS-II requires an input map of ecoregions as biophysical inputs. LANDIS-II ecoregions are continuous or discontinuous areas of the landscape which share climate conditions and soil texture (see section “Model Description”). To generate the map of ecoregions, we used a reference climate dataset and a regional soil map (Table 1). The reference climate dataset reports annual precipitation and mean, minimum and maximum temperature for the period 1971–2000 at a 100 m resolution. Firstly, an unsupervised k-means clustering was applied to the climate dataset to lump together cells with similar climate. MacQueen algorithm was used in this clustering (MacQueen, 1967). We evaluated different numbers of clusters and eventually ended up with four climate regions. The resulting ecoregions agree with our expectations considering topography and the bioclimatic zones found in the study area (Figure 3).

FIGURE 3

Figure 3. Schematic view of workflow followed to produce the ecoregions map. Reference climate data consist of annual precipitation, mean, minimum and maximum annual temperature for the period 1971–2000 at 100 m resolution. Clustering of climate data resulted in four climate regions. Soil data consist of a percentage of silt, sand and clay at 250 m resolution. The intersection of climate regions and soil textures resulted in 27 ecoregions.

The resulting climate regions map was intersected with the soil texture map. The soil map reports the percentage of sand, silt and clay at a 250 m resolution. This map was derived by simply translating the percentages of sand, silt and clay to USDA soil texture categories. The final ecoregions map was therefore produced by overlapping the climate regions and the soil textures maps. This resulted in a total of 28 unique ecoregions defined by both a climate region and a soil type (Figure 3).

PnET-Succession also requires rooting depth for each ecoregion. Ecoregion rooting depth was calculated based on soil depth classes reported in the soil map (0–250 mm, 250–500 mm, 500–1,000 mm, 1,000–1,500 mm, >1,500 mm). The midpoint of each class was used to calculate the most frequent rooting depth for each ecoregion. Precipitation Loss Fraction and Leakage Fraction were given values of 0.6 and 1 for all ecoregions. All other PnET-Succession ecoregions parameters were given default values. A complete dataset containing model inputs is provided (see section “Data availability”).

Calibration of Model Parameters

LANDIS-II and PnET-Succession require a set of parameter values for each species simulated in the landscape. Generic species parameters are required irrespective of the chosen succession extension (e.g., longevity, sexual maturity), while others are required by PnET-Succession (e.g., foliage nitrogen, foliage turnover, minimal, and optimal photosynthetic temperature). We firstly defined an initial set of parameter values from multiple sources. Then, we ran single-cell simulations to verify the species behaviour (e.g., growth, photosynthetic rates) under different conditions. On a single, empty 100 m cell, a single 10-years old cohort of each species was initialized and grown for 200 years preventing establishment of new cohorts. Five replicates were run for each simulation using static monthly averages of temperature, precipitation, PAR and CO₂ (Table 1 and Supplementary Figures 4–6). Baseline climate conditions were used to avoid introducing variability due to fluctuating climate (Gustafson and Miranda, 2019).

The results of these simulations were compared with species biomass estimations for the study area and literature information. Biomass estimations were calculated for P. halepensis, P. nigra, P. pinaster, P. sylvestris, Q. ilex, and Q. pyrenaica based on forest inventory data and allometric equations published in Montero et al. (2005). Simulated versus observed Relative Growth Rates in relation to species biomass were used for comparison as this metric is independent from age, which has a high uncertainty degree in NFI data.

Following calibration guidelines for PnET-Succession (De Bruijn et al., 2014; Gustafson and Miranda, 2019) and evidences from a sensitivity analysis with the same model (McKenzie et al., 2019), we adjusted the most influential species-specific parameter values in an iterative process until the species showed the expected behaviour based on authors’ expertise and observations.

We also evaluated the response of competing assemblages of species typical of the different forest types included in our landscape. These multi-species simulations allow the calibration in relative terms, as well-known competition effects can be assessed and species parameters can be adjusted accordingly. In these simulations, species are established in the cell at the same time and no new establishment is allowed, which is often not the case in natural ecosystems. Thus, the observed species development is the result of their growth traits and their different performance under competition and not due to different establishment strategies or other advantages.

Landscape Simulation

We simulated forest dynamics at the landscape scale incorporating the effect of spatial processes such as dispersal and climate and soil heterogeneity. We initialized the model with the biophysical inputs described above and with the parameters calibrated in the previous step. To verify the influence of shrub communities on simulated forest dynamics, we initialized LANDIS-II with two different vegetation datasets: (1) With shrub communities and (2) Without shrub communities. Since the aim of this study was to initialize the landscape for further experiments with LANDIS-II and PnET-Succession, neither natural nor human-driven disturbances (i.e., fire, harvest) were included in the experiment. We ran five model replicates using baseline climate for 200 years (Table 1 and Supplementary Figures 4–6). We analysed model outputs in terms of temporal patterns of average biomass for each species. Moreover, we mapped and compared total aboveground biomass for selected simulation years across the landscape. All analyses were performed in R version 3.5.3 (R Core Team, 2020) and QGIS 3.10 (QGIS Development Team, 2020).

Results

Spatial Imputation and Initial Vegetation Map

The analysis of 981 NFI plots falling within the study area resulted in eight species having a coverage value higher than 1%: Quercus ilex (25%), Pinus sylvestris (18%), Pinus pinaster (16%), Pinus nigra (16%), Pinus halepensis (16%), Quercus pyrenaica (2%), Populus nigra (1%), and Juniperus oxycedrus (1%). These species were selected to be included in the study. Besides, two extra species—Quercus faginea and Juniperus communis—were also included due to their importance in specific environments. J. communis covers vast areas above the tree line (cushion shape shrubs) (García, 2001) and Q. faginea is also locally abundant.

The analysis of the FM resulted in a classification of active versus inactive cells within the study area. Inactive areas cover 19% of the study area and mainly refer to crops and firebreaks (17.3%). Moreover, active cells were classified as forest and non-forest, which represent 41 and 40% of the study area, respectively.

The intersection of NFI plots and FM polygons defined the species-age assemblage of a total of 789 polygons within the study area. Out of 3,113 polygons with no NFI plot in them, 84% of them were imputed based on equivalent polygons (55% full equivalent, 22% partial equivalent and 7% species equivalent). The remaining 16% of polygons were imputed the species assemblage reported in the FM and the most common age for each of the species.

The use of the VM allowed a more detailed description of active non-forest areas. The analysis of sparse trees in non-forest polygons increased by 4% the surface of the study area where tree species are present. Even though it may seem as a small portion of the landscape, these sparse trees can represent important seed sources when long-term forest dynamics are simulated. Moreover, the VM analysis allowed the inclusion of shrub communities in the landscape’s cells, which can affect the shade conditions as well as water availability.

Figure 4 shows the initial vegetation map as a result of combining the presence of forest trees, sparse trees and shrubs communities. Several portions of the landscape are covered by shrublands and Juniperus spp., surrounded by a mosaic of Pinus spp., Quercus spp. and mixed Pinus-Quercus forests, with a minor presence of Populus nigra.

FIGURE 4

Figure 4. Map of initial vegetation conditions including shrub communities. Each category represents a community dominated by one or more tree species, where shrubs and Juniperus spp. may also be present.

Calibration of Model Parameters by Means of Site-Level Simulations

Estimating species parameter values was a complex task due to the amount of information sources required to cover all of them. Table 2 reports the most important parameters required by LANDIS-II and PnET-Succession. A detailed explanation of all parameter values, sources and rationales for their adoption is provided in Supplementary Table 2.

TABLE 2

Table 2. Species parameter values for LANDIS-II and PnET-Succession.

Biomass estimations for P. halepensis, P. nigra, P. pinaster, P. sylvestris, Q. ilex, and Q. pyrenaica were compared with results from single-cell simulations initialized with individual species. Results of biomass and photosynthetic rates of species simulated individually on a single cell are reported in the Supplementary Figures 1–3. Simulated results are within the range of estimations, although estimations are highly variable among plots.

Among the more than 70 of single-cell simulations that were performed, we chose to show here two representing typical species assemblages of low and high altitude forest types (Figure 5). In the low altitude forest type, composed of two pine and two oak species, our results indicated a dominance of Q. ilex over Q. faginea, P. halepensis, and P. pinaster (Figure 5A). Cohorts of the two latter species were simulated to die by year 140 and 190, respectively, since they approached their longevity (pink and blue lines). Q. ilex clearly dominated Q. faginea but it did not fully outcompete it. The advantage of Q. ilex in this forest community seems to be related to its capacity to start photosynthesizing earlier in spring than the other species (Figure 5B).

FIGURE 5

Figure 5. Mean simulated annual total biomass (wood, roots and foliage) (A,C) and monthly net photosynthesis (B,D) for the two chosen forest types. Shaded area refers to standard deviation. Values of net photosynthesis are averages between years 50–75, with error bars indicating standard deviation. Note that lines connecting the points serve for illustrative purposes only, since PnET-Succession simulates photosynthesis at monthly time scale (not daily).

In the high altitude forest type, P. sylvestris, P. nigra, and Q. pyrenaica coexisted along the simulation, although P. sylvestris built higher biomass compared to the other two species (Figure 5C). The advantage of P. sylvestris was related to its higher photosynthetic rate from the beginning of the season, while Q. pyrenaica, being a deciduous species, increased its photosynthesis more gradually after having built foliage biomass (Figure 5D). Generally, we found that PnET-Succession reproduced reasonably well the bimodal growth patterns of Mediterranean tree species, mostly occurring during spring and fall instead of summer which is characterized by a lack of precipitations.

Landscape-Level Simulation

Both simulations with and without shrubs showed a trend to increase the average biomass of all tree species during the first years of the simulation and stabilization afterwards (Figure 6). In the simulation including shrubs, this increase was slower than in the simulation without shrubs, with a faster increase in the first 50 years. In both simulations the total average biomass at stabilization was around 100 tn ha^–1.

FIGURE 6

Figure 6. Average biomass of species through time for the whole landscape. Note that Junipers (J. communis and J. oxycedrus) and shrub communities (short, medium, and tall) have been grouped together.

We found differences among species in terms of simulated biomass growth. Among pine species, P. sylvestris and P. nigra were those with a higher increase of average biomass in both simulations. Biomass of P. halepensis and P. pinaster increased during the first 50 years, followed by a decline and disappearance from the landscape toward the last decades of the simulation. Among oaks, biomass of Q. ilex increased notably, while Q. pyrenaica increased at lower rates and stabilized after about 100 years. Biomass of Q. faginea and Pop. nigra had similar trends, increasing slightly in the early years and declining afterwards, but still maintaining presence at low levels of biomass. Junipers slightly increased their biomass during the first 30 years and then strongly declined. These species show similar patterns in both simulations, but under the simulation without shrubs we observe a steeper increase of biomass during the first years. Shrub communities increased their average biomass during the first 40 years and declined afterwards.

Initial total landscape biomass (time 0) was similar in both simulations since shrub communities accounted for low levels of biomass (5.1% of total biomass). At time 50, the area with high biomass is wider in the simulation without shrubs (Figure 7). This pattern was observed at time 100 too, although the difference between the two simulations was smaller (see total biomass at each time step in Figure 7). At time 150 and 200, both simulations showed a similar quantity and distribution of the biomass across the landscape, with slightly higher values of biomass in the simulation without shrubs than in the one with shrubs. The biomass distribution within the study area followed the altitudinal gradient, with higher biomass found at medium elevations, specially at time step 150. At the beginning of the simulation, areas at high elevation (above c.a 2,000 m a.s.l.) show the smallest values of biomass relatively to the rest of the study area. This pattern remained by the end of the simulation in both cases.

FIGURE 7

Figure 7. Total aboveground biomass (tn ha-1) at years 0, 50, 100, 150, and 200 for the two simulations, as indicated in the scale bar on the right. Number in each panel indicates the total aboveground biomass (Tg) for the whole landscape at each time step.

Discussion

We provide detailed step-by-step example to initialize, calibrate and set up a forest landscape model. Our work could help other potential users to better understand what is required to start applying such models. Thus, our fully documented methodological process represents a step forward toward the transparent application of forest landscape models in regions without prior application. We also made available a high-resolution map of vegetation conditions and calibration details for a large mountainous landscape in the European Mediterranean area, together with the input data and scripts used in the process. Our landscape level simulations reveal distinct dynamics among species according to their competitive potential and simulated intra-annual growth. These results also indicate that shrub communities shall be considered in forest landscape models as they have the potential to influence forest dynamics by delaying growth and expansion of tree species in Mediterranean ecosystems.

Spatial Imputation and Initial Vegetation Conditions

The selection of species to be included in modelling studies is typically done by analysing data from terrestrial plot measurements such as forest inventories (e.g., Wang et al., 2019). However, for certain ecosystems, limiting this inspection to forest inventories has some disadvantages. Forest inventory data commonly report information on tree species, neglecting important functional groups (in our case, J. communis and shrubs communities). Moreover, selection is usually based on variables such as basal area or stand coverage. This selection may result in the exclusion of species that may not be abundant at the landscape level but whose presence is crucial at the local scale. For example, Q. faginea would have been initially excluded from our study since its coverage falls below the 1% of the study area (0.03%). However, this species is found at high abundance in some locations and has a higher susceptibility to summer drought than Q. ilex. Therefore, the dynamics and distribution of Q. faginea are expected to be highly affected by climate change (Quero et al., 2006). Similarly, J. communis was included due to its importance in areas above the treeline. In these areas, J. communis is susceptible to interact with tree species under climate change conditions by limiting or facilitating uphill migration of tree species. Therefore, even though forest inventory data are important resources for generating inputs for FLMs, we recommend combining such datasets with other information sources (e.g., vegetation maps, local studies) according to the forest ecosystem under consideration and study focus—in our case, climate change applications—. Considering factors such as relevance in specific habitats or cultural importance (e.g., Lucash et al., 2019) may result in inclusion or exclusion of some species.

Our method, which combined multiple information sources allows the inclusion of fine-coarse information such as the presence of sparse trees and shrub communities. The inclusion of sparse trees in the initial vegetation map is important in the Mediterranean mountains, since we simulate succession in pine plantations, artificially created stands where regeneration is highly affected by seed dispersal from adjacent patches of native vegetation (González-Moreno et al., 2011; Navarro-González et al., 2013). Although shrubs’ biomass is generally low, they can shade the forest floor and therefore influence simulated light and soil water dynamics, as well as affect establishment (further discussion below). Moreover, shrubs play a key role in other processes such as fire dynamics due to their role as fuel in forest fires (Syphard et al., 2006).

In this study, we have used a categorical methodology for the spatial imputation, contrary to commonly used methods based on geographic or data space distance (e.g., Wilson et al., 2012; Ohmann et al., 2014; Duveneck et al., 2015). Our study area is a highly anthropized mosaic of different land uses and management regimes, and therefore imputation based solely on distance might not have been an appropriate criterion, as indicated by Duveneck et al. (2015). For example, the most common forest type in our landscape, pine plantations, are the result of past forestry policies which were applied almost simultaneously all over the region. Therefore, a high similarity is expected between stands, regardless of the physical distance between them. This is a common situation in other Mediterranean ecosystems (European Environmental Agency [EEA], 2016), and therefore this approach could also be used in such cases.

The uncertainties of this methodology are inherently related to the uncertainties of the information sources. Besides, an additional source of uncertainty in the initial vegetation conditions map is related to the assignation of age to cohorts. Tree age data is rarely available at single-tree level in forest datasets, as reliable tree age estimations are resource-consuming and often invasive (Fazan et al., 2012; Rohner et al., 2013). Thus, modellers are commonly forced to assign each tree (or cohort) to an age class inferring it from available measurements (e.g., diameter, height, average stand age; Abrams, 1985; Rozas, 2003). In the case of LANDIS-II, a cohort’s age determines simulated biomass, which in turn influences light and water availability at cell-level. Therefore, assumptions in the age assignation process may not be so relevant as long as the relative difference between species reflects real conditions. Considering this, we created correspondence rules between tree age and diameter and height for each of the simulated species and age classes. Thus, the growth pattern assigned to each of the species is relative to the other species, reflecting the differential access to water and light by each species-age. This methodology could be further improved by considering additional local yield tables and observations. Nevertheless, by documenting such correspondence rules the model inputs generation is significantly more transparent and reproducible than other LANDIS-II studies (e.g., Mina et al., 2021).

Our goal was to increase the reproducibility of model input generation by ensuring high transparency and detail in the process description. The methodology presented in this study does not necessarily introduce new methods compared to previous studies (e.g., He et al., 1998) but rather it highlights all aspects of the process, which we believe could be of great benefit for beginner modellers to set up applications in new landscapes. Firstly, considering multiple information sources at plot- and polygon-level may be necessary. In our case, multiple vegetation maps were required to consider forested areas, dispersed trees and shrub communities. Secondly, the collected information likely requires processing and transformation, which may introduce assumptions (e.g., age assignation). Thirdly, the selection of the appropriate spatial imputation method should consider the study area characteristics (e.g., coetaneous patches of pine plantations) and available information.

Calibration of Model Parameters and Site-Level Simulations

Calibration of model parameters was performed by running simulations and testing long-term species dynamics and competitive interactions at site-level. The obtained results were iteratively assessed to adjust parameter values until the species showed their expected growing patterns. Experiments at site scale using landscape models have also been used before to analyse the influence of different factors on model outputs by avoiding the high complexity resulting from large landscapes simulations (Gustafson et al., 2017, 2018).

During the calibration phase, species parameters were adjusted to ensure that the model simulates realistic species biomass estimations. Species biomass data derived from field observations (e.g., growth-and-yield sites, inventories and permanent plots) are usually highly variable as they differ depending on multiple factors (e.g., location, stand development, site index). Moreover, biomass values often have a high degree of uncertainty, since they are commonly estimated based on general allometric equations from other measured variables (e.g., diameter, height, wood density) (Forrester et al., 2017). Therefore, such comparisons should be interpreted with caution. In this study we used such estimations to ensure that the simulated biomass falls within realistic ranges rather than adjusting parameter values to match the exact values (see Supplementary Figure 3). With this approach, we calibrated the most relevant species for this study (P. halepensis, P. nigra, P. pinaster, P. sylvestris, Q. ilex, and Q. pyrenaica).

Besides, we review here other studies which provide biomass estimations for some of our species. The dynamics of Q. pyrenaica stands have been studied by Santa Regina (2000), who estimated its biomass in four plots in north-western Spain. Our simulations have a high degree of agreement for foliage biomass, while certain overestimation remains for wood biomass. This difference, nevertheless, can be justified as in the plots studied by Santa Regina (2000) the presence of shrubs could be reducing Q. pyrenaica productivity. Río and Sterba (2009) studied the productivity of mixed stands of P. sylvestris and Q. pyrenaica. They found that although P. sylvestris is less productive when growing in mixed stands, the reduction in productivity is smaller than the reduction in occupied area. Accordingly, our simulations show a decrease in productivity of P. sylvestris when growing together with other species such as Q. pyrenaica, but it remains as a highly productive species.

Other methodologies have been applied for model calibration (e.g., Duveneck et al., 2017; Cassell et al., 2019; Mina et al., 2021). As an example, Duveneck et al. (2017) used data from flux towers within New England (United States) to calibrate PnET parameters. The application of some methodologies over others usually responds to the availability of data for the simulated area. In this sense, the lack of available biomass accumulation curves limited the application of more exhaustive calibration methodologies. Our calibration could therefore be improved if additional data sets become available, such as high-resolution biomass measurements or growth rates based on flux towers measurements.

Our results clearly show that PnET-Succession reproduces the characteristic bimodal growth observed in Mediterranean species (Larcher, 2000). This growth pattern is the result of dry summer conditions, which impose a limitation for growth to numerous species. Thus, Mediterranean species often show two peaks of productivity through the year, in late spring and fall (Camarero et al., 2010; Gutiérrez et al., 2011). Q. ilex bimodal growth has been studied in detail by Gutiérrez et al. (2011). They report asymmetrical radial increment peaks in May and September for a coastal location in north-eastern Spain, with high plasticity dependent on climatic conditions and most of the growth occurring during the first growing phase. The simulated growth of Q. ilex reproduces this pattern both when the species is simulated growing alone or in association with other species. Modelling such growth interannual variability remains a challenge for forest models, for which improvements are proposed to include such processes (Mina et al., 2016). In the case of PnET-Succession, the high temporal resolution (monthly scale) and mechanistic approach allows reproducing such patterns and therefore makes this model highly suitable for applications in Mediterranean systems.

Distinct dynamics were observed according to the species competitive potential and simulated intra-annual growth. Species leaf habits is one of the factors influencing species competitiveness. Among the simulated species considered in this study, only Q. faginea, Q. pyrenaica, and Pop. nigra are deciduous while all other species are evergreen. Deciduous species have a higher potential for net photosynthesis in PnET-Succession (De Bruijn et al., 2014) but they are generally less shade and drought tolerant and need to spend more energy in building foliage biomass than evergreen species. The trade-off between benefits and losses caused by different leaf habits explains the coexistence of species with different strategies (Escudero et al., 2017), as it is clearly observed in the growth patterns of Q. ilex and Q. faginea.

The calibration of non-tree species—junipers and shrubs—was particularly challenging. Although we believe including both junipers (J. communis and J. oxycedrus) in the initial vegetation map was important in our target ecosystem, the lack of information and reliable data on these species limited a finer calibration. As a result, both species were assigned with very similar parameter values, their simulated behaviour was almost identical and thus they were grouped together. In the case of shrub communities, since they represent functional groups rather than single species, calibration was achieved mostly by extrapolation, comparison with similar studies (e.g., Cassell et al., 2019) and tuning according to expected simulated behaviour relatively to tree species. Their ecological role, for the sake of this study, was mainly as shade providers and competitors for establishment, and therefore our main objective in this sense has been oriented to ensure survival and growth beneath the tree canopy. For further applications where the role of non-tree species is more relevant (e.g., fire dynamics, facilitation), their parameterization shall be improved to better reflect differences between shrubs species or functional types.

Landscape-Level Simulation

Our results show that the shrub communities influence the forest dynamics by delaying the growth and expansion of tree species. In our simulations, we found that shrubs reduce tree species establishment. However, shrubs are known to serve as nurse plants, favouring tree seedling growth by amelioration of adverse dry conditions and protection against herbivory (Castro et al., 2004; Gómez-Aparicio et al., 2005; Prévosto et al., 2020).

Modelling the role of understory in forest succession has also been investigated by Thrippleton et al. (2016) with the LandClim model. Similarly to our results, they found delayed establishment of trees when herbaceous understory was abundant. Moreover, both Thrippleton et al. (2016) and our results show that shrubs are unable to establish under dense canopies, eventually declining and even disappearing from the landscape. This agrees with field observations: weak regeneration of shrubs and tree species under dense canopies biomass even when the sapling bank is present (Mendoza Sagrera, 2008). Also, the decline of shrubs in our simulations was likely not realistic, since small-scale perturbations creating patches where light availability increases and shrub communities thrive (e.g., due to fire, Leverkus et al., 2014) were not included in our experiment. However, the interactions occurring between trees and non-trees species and abiotic factors such as fire justify the need to include these communities in these kinds of applications (e.g., Loudermilk et al., 2013).

The increase of forest biomass observed in our simulations at landscape level was somewhat expected, since we did not include harvesting or natural perturbations (e.g., fire, pest outbreaks). Additionally, simulations were run with the same climate inputs used for calibration (baseline climate), thus potential impacts of changing climate (e.g., higher temperatures, extended drought, CO₂ fertilization) were not considered. At medium altitudes, where most forest stands are found, the growing tendency in biomass was likely a result of pine plantations being relatively young at the beginning of our simulations. Pine plantations in Sierra Nevada and Sierra de Baza have been showing in latest years signs of decay as a result of increasing drought stress and intense interspecific competition (e.g., Sánchez-Salguero et al., 2010, 2012a, b). These mortality events could not be observed in our study since increasing drought was not considered in our climate inputs.

Even though there was a general growing trend in biomass, the dynamics among species differed. P. halepensis and P. pinaster biomass declined as they approached longevity and seemed unable to regenerate beneath the canopy, while Q. ilex kept growing. P. sylvestris and P. nigra, however, were able to coexist with Q. pyrenaica, and therefore they remained present in the landscape and even increased their biomass through time. Factors influencing the establishment of pines and oaks in Sierra Nevada and Sierra de Baza have been profusely studied (e.g., Mendoza et al., 2009a, b; Gómez-Aparicio et al., 2009; Herrero Méndez, 2012), showing limitations to recruitment mainly due to high post-dispersal seed predation rates and dry summer conditions. Our simulations also reflect poor establishment of some species, probably related to dry summer conditions and species shade tolerances being higher than in reality. In the case of junipers, they are probably limited by a short growing season due to their minimum and optimal photosynthesis temperatures and low biomass levels which prevents them from accessing light and water.

Modelling Aspects, Limitations and Future Research

The selection of the model used in this study was based on the flexibility and potentiality of this model. LANDIS-II is a well-documented model supported by an active scientific community, which ensures help and support for new modellers. Moreover, several model extensions are available for a wide variety of future applications (e.g., harvest, fire, wind, biological agents). To simulate forest dynamics, we specifically chose PnET-Succession as it simulates ecological succession with a more mechanistic approach than past extensions. As mentioned above, by simulating growth at monthly resolution, intra-annual growth variability can be properly captured, which is of crucial importance for simulations under climate change and applications in the Mediterranean area.

However, the model is still limited in some aspects which could be further improved. First, the need to provide cohort age certainly adds uncertainty to the initial conditions. Given that PnET-Succession uses age to determine initial cohorts’ biomass, other variables such as diameter, basal area or height—often available from field observations—could be used instead, reducing uncertainty. Even though latest PnET-Succession version allows providing initial biomass, this variable often needs to be estimated with allometric equations. Besides, since cohort’s biomass is highly dependent on tree density, tree density could also be incorporated into the model in some way to improve model representation of competition effects (but see Wang et al., 2014). This would be especially relevant in ecosystems such as pine plantations, where extremely high tree densities cause inter- and intraspecific competition to be an important factor explaining vegetation dynamics.

The grouping of cells into homogeneous climatic conditions (ecoregions) prevents a fine-coarse capturing of topographic influence on climate, which is a relevant issue in mountainous areas. The definition of such ecoregions is commonly done by clustering average climate information to define homogeneous climate regions. These homogeneous regions are assigned with climate time series typically obtained as the average among all cells within the ecoregion. By doing so, the influence of topography on climate is not well captured and extreme conditions, such as the ones found in the mountain peaks, are smoothed down. In our case, the increase in biomass at high altitudes was probably due to the fact that temperature was not limiting enough the establishment of species beyond their natural treeline. In a large model comparison study, Petter et al. (2020) also found that no clear treeline emerged from LANDIS-II simulations, in which the area above vegetation was binned into a single ecoregion, compared to other models in which each cell has its own climatic condition. The use of continuous maps for climate conditions (i.e., each cell assigned its own environmental and climate condition) could improve how PnET-Succession simulates the effect of topography on climate and thus on vegetation changes. Other forest landscape models developed for mountain environments already make use of continuous maps explicitly accounting for topography (e.g., LandClim; Schumacher et al. (2004)). This feature is currently under development for LANDIS-II, which was initially designed to efficiently simulate large landscapes (more than 1 million interacting cells).

Even considering the above-mentioned limitations, the model showed a great potential for a wide range of applications in the Mediterranean area. The transparent initialization of the model and the documented calibration can serve as a guide for new users, encouraging the application of forest landscape models. Besides, testing different initialization data has allowed us to confirm the importance of shrub communities in the forest dynamics within the study region. Further experiments will analyse the vegetation dynamics under natural perturbations such as fire or defoliators. Moreover, the inclusion of climate change and silviculture will allow us to explore future forest dynamics, and, by doing so, deliver management recommendations to promote ecosystems adaptation to global change.

Data Availability Statement

The datasets and scripts generated for this study can be found in https://doi.org/10.5281/zenodo.4584475.

Author Contributions

MS-M: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, software, visualization, writing – original draft, writing, review, and editing. MM: conceptualization, methodology, supervision, validation, writing, review, and editing. PS: formal analysis, writing, review, and editing. RN-C and JQ: validation, writing, review, and editing. FB-G: conceptualization, resources, writing, review, editing, supervision, methodology, and funding acquisition. All authors contributed to the article and approved the submitted version.

Funding

MS-M was supported by the Spanish Ministry of Science and Innovation (FPU predoctoral grant and the project PROPIFEN PGC2018-101773-B-I00). MM acknowledges a grant from the Swiss National Science Foundation (project n.175101) and funding from the Canada Research Chairs Program. PS was funded by the project “Ecología Funcional de los Bosques Andaluces y Predicciones Sobre Sus Cambios Futuros” (For-Change) (UCO-27943) by Junta de Andalucía (Spain), the project “Funcionalidad y servicios ecosistémicos de los bosques andaluces y normarroquíes: relaciones con la diversidad vegetal y edáfica ante el cambio climático” by “Ayudas a la I+D del Plan Andaluz de Investigación, Desarrollo e Innovación (PAIDI) 2020,” Junta de Andalucía (Spain) and European FEDER funds. JQ and RN-C were funded by the project ESPECTRAMED (CGL2018-86161-R) from Spanish Research Agency, Ministry of Science and Innovation. RN-C was funded by the projects ISOPINE (UCO-1265298) and SilvAdpat Network RED2018-102719.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We acknowledge the E-OBS dataset from the EU-FP6 project UERRA (http://www.uerra.eu) and the Copernicus Climate Change Service, and the data providers in the ECA&D project (https://www.ecad.eu). We thank Robert Scheller, Brian Miranda, Rafael Villar, and Núria Aquilué for useful insights on project conceptualization during the initial steps of this work. We also thank two reviewers for the comments and suggestions made to earlier versions of this manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2021.653393/full#supplementary-material

Footnotes

^ www.landis-ii.org/

References

Aber, J. D., and Federer, C. A. (1992). A generalized, lumped-parameter model of photosynthesis, evapotranspiration and net primary production in temperate and boreal forest ecosystems. Oecologia 92, 463–474. doi: 10.1007/bf00317837

PubMed Abstract | CrossRef Full Text | Google Scholar

Abrams, M. D. (1985). Age-diameter relationships of quercus species in relation to edaphic factors in gallery forests in Northeast Kansas. For. Ecol. Manag. 13, 181–193. doi: 10.1016/0378-1127(85)90033-7