Structural uncertainty and uncertainty management in four common Land Use Cover Change (LUCC) model software packages. A comparison

Research on the uncertainty of Land Use Cover Change (LUCC) models is still limited. Through this paper, we aim to globally characterize the structural uncertainty of four common software packages (CA_Markov, Dinamica EGO, Land Change Modeler, Metronamica) and analyse the options that they offer for uncertainty management. The models have been compared qualitatively, based on their structures and tools, and quantitatively, through a study case for the city of Cape Town. Results proved how each model conceptualised the modelled system in a different way, which led to different outputs. Statistical or automatic approaches did not prove to be more successful that userdriven approaches. The available options for uncertainty management vary depending on the model. No model pays enough attention to the communication of uncertainties.

▪ Model operation uncertainty. The uncertainty that arises from the accumulation and interaction of uncertainties propagated through the model. These sources of uncertainty in LUCC modelling exercises have been addressed widely in literature. Many papers assess the output uncertainty of specific exercises (Memarian et al. 2012;Ligmann-Zielinska 2013), whereas others focus on specific sources of uncertainty, like parameter uncertainty (Dietzel and Clarke 2004;Conway 2009;García et al. 2011;Van Vliet et al. 2013b;Houet et al. 2015;Confalonieri et al. 2016;Grinblat et al. 2016;Liao et al. 2016) and process variability uncertainty (Kok and Van Delden 2009;Verburg et al. 2013;Alexander et al. 2015;Maier et al. 2016).
Less common is the research about epistemological uncertainties of LUCC models and, specifically, about the model's structural uncertainty. Some studies address specific topics related to this issue, like the different procedures to calculate change potential and change allocation (Riveira and Maseda 2006;Lin et al. 2011;Pérez-Vega et al. 2012;Camacho Olmedo et al. 2013;Shafizadeh-Moghadam et al. 2015). Ferchichi et al. (2017b) propose a framework to quantify structural uncertainty of LUCC models based on probabilistic theory and sensitivity analysis. However, there is a lack of studies that This manuscript has been submitted for publication in "Environmental Modelling & Software". Subsequent versions of this manuscript, after acceptance and review by the journal, may have slightly different content. If accepted, the final version of this manuscript will be available via the `Peer-reviewed Publication DOI' characterize the overall structural uncertainty of available LUCC model software packages.
Today, there is a large availability of standard model software packages to simulate different spatial dynamics ) and their use is ever-increasing (Wickramasuriya et al. 2009;Chaudhuri and Clarke 2013;Leija et al. 2021). Information about the uncertainty associated to the use of these software packages is not usually widespread and no paper analysing their structural uncertainty has been found in the literature. However, knowledge about the uncertainty of these models is essential to engage planning agents and spread their use in real-world problems solving (Yeh and Li 2006;Batisani and Yarnal 2009;Sohl et al. 2016), improving the understanding and characterization of these tools.
Through this paper, we aim to fill the previous research gap. With that objective in mind, we compare four standard LUCC model software packages and characterize their structural uncertainty. Model comparison has been proposed by several authors as a way to assess the structural uncertainty (Kelly et al. 2013;Uusitalo et al. 2015) and has been usually employed as an useful approach to better characterize and understand the available software packages (Toro Balbotín 2014; Mas et al. 2014;Aguejdad et al. 2016;. Through the comparison, users may be aware of the pros and cons a specific approach may entail and which may better fit their interests.

Model software packages
We compared four standard LUCC model software packages: CA_Markov (Eastman and Toledano 2018a), Dinamica EGO (Soares-Filho et al. 2002, 2009Rodrigues and Soares-Filho 2018), Land Change Modeler (Eastman 2015;Eastman and Toledano 2018b) and Metronamica (RIKS 2012;Van Delden and Vanhout 2018). They are a good representation of LUCC modelling software commonly employed in literature (Santé et al. 2010;Kamusoko 2012;Eastman and Toledano 2018b;Ferreira et al. 2019). A brief description and graphic representation of each one is provided in the Annex 1. Annex 2 includes a comparative table to evaluate the differences among models. This manuscript has been submitted for publication in "Environmental Modelling & Software". Subsequent versions of this manuscript, after acceptance and review by the journal, may have slightly different content. If accepted, the final version of this manuscript will be available via the `Peer-reviewed Publication DOI'

Model comparison
Model structures were compared according to the following aspects: change potential calculation, quantity of changes estimation, allocation of changes, pattern simulation, and validation and outputs. The comparison followed a both qualitative and quantitative analysis of each software (Fig. 1).

Figure 1. Conceptual chart of the followed methodology for models' comparison
This manuscript has been submitted for publication in "Environmental Modelling & Software". Subsequent versions of this manuscript, after acceptance and review by the journal, may have slightly different content. If accepted, the final version of this manuscript will be available via the `Peer-reviewed Publication DOI' Through the first approach (2.2.1) we compared the system conceptualization adopted in each model and the available options they offer for uncertainty management. Through the second approach (2.2.2), we compared model outputs for the same case study to assess the differences that come from the model structure. This case study, which is explained in detail in annex 2, consisted of a model set up with a spatial resolution of 100m grid cells for the city of Cape Town and the period 1990-2013.

Qualitative model characterization
We reviewed the main characteristics of the different software packages. This information was obtained through the model's documentation and based on the deep experience of authors with the analysed software. This allowed to understand the conceptualization of each model and the different tools they provide for uncertainty management.

Model outputs assessment
The four models were calibrated for the same case study following the approach described in annex 3. Outputs generated by each model were then analysed and compared. This allowed to analyse which differences between simulations came from the use of different model structures.
For the comparison of output maps, we differentiated between soft-classified and hardclassified maps (Camacho Olmedo et al. 2013). The first ones, which we will also refer to as change potential maps, show the probabilities of change to a specific category. Hardclassified maps, which are the final land use maps simulated by the models and we refer to as simulation results, assign every pixel to a specific category and, therefore, show states instead of probabilities.
Agreement between change potential maps obtained through different methods of change potential calculation was measured through the Pearson correlation coefficient (Mukaka 2012). To this end, transition potential maps to the same category in Dinamica and LCM were aggregated and compared to the land use potential maps for that category in CA_Markov and Metronamica. These last maps were masked to only account in the comparison for those areas where the transitions considered in Dinamica and LCM could take place.
Simulation results from different models were compared by means of standard crosstabulation techniques, Kappa Simulation (Ksim) and a set of spatial metrics calculated at This manuscript has been submitted for publication in "Environmental Modelling & Software". Subsequent versions of this manuscript, after acceptance and review by the journal, may have slightly different content. If accepted, the final version of this manuscript will be available via the `Peer-reviewed Publication DOI' the class level: number of patches, patch mean size and standard deviation and proportion of like adjacencies. Ksim evaluate the agreement between the changes simulated by each model compared to the agreement that is expected by chance (Van Vliet et al. 2011).
Spatial metrics characterize the shape and of patches and the way they are allocated on the map, that is, the maps´ patterns. A patch is a group of neighbour pixels with the same value (Botequilha Leitao et al. 2006). The proportion of like adjacencies inform about the aggregation of cohesion between patches of the same class (Mcgarigal 2018). That is, how aggregated or fragmented are the patches that make up a class.
Each model was executed 20 times under the same parameters and conditions to assess the intra-model output variability. Only outputs from the first executions were employed for the previous assessments, whereas outputs from the remaining 19 executions were only employed at this step to assess the agreement between model executions. Agreement between change potential maps obtained through the same production method for each model was assessed by calculating the average standard deviation of the pixel values across the 20 outputs. Agreement between simulations was assessed through KSim and cross-tabulation, as in the inter-model comparison.

LUC and LUCC explanatory factors
Each model imposes limitations on the number and nature of the LUC or LUCC explanatory factors to be considered. Metronamica only comprises four factors: neighbourhood interactions, accessibility, suitability and zoning. Notwithstanding, the user can choose how many base maps he or she likes to include in the different factors and even rule out some of them by modifying the transition potential formula that guides the change potential map creation.
CA_Markov, Dinamica EGO and LCM do not impose any limits on the number and thematic of the factors to be considered. It is the user who should decide this based on his knowledge and information about the modelled system. In addition, LCM includes specific tools to transform the model's factors, associating them with specific transitions of change of land uses / covers, to improve the statistical explanatory power of the drivers.
When using this option, as far as the model driver is transformed based on LUC data, the This manuscript has been submitted for publication in "Environmental Modelling & Software". Subsequent versions of this manuscript, after acceptance and review by the journal, may have slightly different content. If accepted, the final version of this manuscript will be available via the `Peer-reviewed Publication DOI' model builds on the assumption of temporal stationarity, i.e. the continuation of past trends to the future, which may not be real.
Metronamica allows to dynamically change the explanatory factors throughout the temporal extent of the simulation. Dinamica EGO and LCM can also work with dynamic LUCC explanatory factors. This allows to manage the uncertainty to which static approaches do not give an answer. CA_Markov does not offer any option to this end.
In our study case, all models were calibrated according to a similar set of factors: accessibility to roads, train stations and economic hubs; slope; median price values; protected areas; and land use zoning. LCM and Metronamica also included factors that accounted for the power of attraction of simulated categories to attract new land uses.
These proved to be the most important factors of their respective applications. In CA_Markov and Dinamica EGO these factors were not considered, as the CA component of these models already accounts for that attraction rules in the allocation of changes step (see section 3.4).

Two different types of change potential maps
Two types of change potential maps can be distinguished: Transition Potential (TP) maps LUP maps indicate the preference of each land use class to occupy any location of a study area based on a set of drivers defined by the user. TP maps indicate the potential of a set of defined land uses to transition to another set of land uses. Models that use LUP maps are more flexible than models that use TP maps. Although transitions can take place from any class to every other class on TP maps, they are usually restricted to the more meaningful. In LUP maps, changes can take place anywhere in the study area and are not limited to a specific set of transitions.
The change potential calculation of TP models usually relates historical changes and driving forces. In a similar vein, LUP maps are usually created based on an understanding of historical changes. According to Camacho Olmedo et al. (2013), when creating these maps, we implicitly account for all past changes, including remote past. However, in the case of LUP maps the user may decide to do not take that historical information into This manuscript has been submitted for publication in "Environmental Modelling & Software". Subsequent versions of this manuscript, after acceptance and review by the journal, may have slightly different content. If accepted, the final version of this manuscript will be available via the `Peer-reviewed Publication DOI' account. That is not always the case of TP maps. Dinamica EGO and LCM require LUC maps at 2 different time points to calculate the TP maps. Although Dinamica EGO allows for manual adjustment of the statistical relation found by the model through the WoE method, LCM does not offer any option to this end.

Change potential calculation
Each model calculates the change potential though a different procedure (Fig. 2).
Nonetheless, the four models offer the user some flexibility regarding the change potential calculation to account for the model structure uncertainty that this step encompasses. This manuscript has been submitted for publication in "Environmental Modelling & Software". Subsequent versions of this manuscript, after acceptance and review by the journal, may have slightly different content. If accepted, the final version of this manuscript will be available via the `Peer-reviewed Publication DOI' methods trust automatic or statistical procedures as defined by Van Vliet et al. (2016) to find out the relation between changes and drivers of change. To this end, they use a sample of pixels as training and then, in the case of Neural Networks and SimWeight, the inferred relations are compared to a set of validation pixels. As far as the analysis sample varies with each model run, the inferred relations change with the sample as well (Kim 2010), although this variation is low ( Table 1). The logistic regression procedure allows the user to employ all the pixels in the analysis and, therefore, avoid this possible uncertainty.
Dinamica EGO makes use of the Weights of Evidence (WoE) to calculate the change potential maps. Soares-Filho et al. (2013) developed a Genetic Algorithm that allows the user to refine the change potential calculated through the WoE method. The software also allows the user to manually edit the obtained weights. It can even use external maps produced through other methods to bypass the incorporated methods.
The manual adjustment of the weights automatically obtained allows the user to account for some of the uncertainties that the data, calibrations periods, etc can entail. However, this adjustment may have a great impact on the obtained maps. In our modelling exercise, change potential maps obtained with automatic and adjusted weights showed big differences (Table 1).
CA_Markov does not integrate a specific method for change potential calculation, although the model help advises to employ the Multicriteria Evaluation (MCE) implemented in TerrSet as the standard tool for this purpose. Accordingly, the uncertainty of the model's change potential maps will be related to the chosen method and its specificities. If choosing MCE, which relies on user or expert knowledge, it will be very dependent on the uncertainty of that knowledge and the criteria included. Once the maps are obtained through the selected method, CA_Markov does not make any modification to them.
In Metronamica, the change potential map is calculated through a formula that combines a series of input data manually adjusted by the user. The user can also edit the formula, which gives him the chance to account for the model structure uncertainty. However, he can also introduce new sources of uncertainty by doing so. Some variability between change potential maps is found at each model execution because of the random component that is included in the creation of this map in Metronamica (Table 1), conceptualized as a factor that accounts for the uncertainty of human behaviour. It is This manuscript has been submitted for publication in "Environmental Modelling & Software". Subsequent versions of this manuscript, after acceptance and review by the journal, may have slightly different content. If accepted, the final version of this manuscript will be available via the `Peer-reviewed Publication DOI' possible to run the model in a deterministic way as well, although the developers do not recommend this due to the inherent uncertainty in land use dynamics.
The correlation between change potential maps produced from different models is usually low, only higher in the case of the maps produced through procedures that involve user knowledge and intervention (Table 1). For example, the correlation between maps produced by the manual adjustment of the Dinamica's Weights is high with maps produced by CA_Markov and Metronamica. The correlation between maps obtained through different automatic or statistical methods, if not using the same software, is very low.
This manuscript has been submitted for publication in "Environmental Modelling & Software". Subsequent versions of this manuscript, after acceptance and review by the journal, may have slightly different content. If accepted, the final version of this manuscript will be available via the `Peer-reviewed Publication DOI' Table 1. In diagonal and grey, average standard deviation of change potential maps produced through the same method after 20 model executions for the transitions to residential areas. Off-diagonal cells show the Pearson correlation coefficient between change potential maps produced through different methods for the transitions to residential areas. CAM: CA_Markov; Metro 0: Metronamica with random factor = 0; Metro 0.5: Metronamica with random factor = 0.5; LG SY10: Logistic Regression in LCM with a Systematic Sampling = 10%; LG SY100: Logistic Regression in LCM with a Systematic Sampling = 100%; LG ST10: Logistic Regression in LCM with a Stratified Sampling = 10%; LG ST100: Logistic  points. In these cases, the uncertainty of input maps will be transferred to the obtained probabilities. Nonetheless, the three models allow the user to manually modify the obtained Markov probabilities from input data, and, therefore, to account for some of the uncertainties associated to input data if wanted.
The Markov probabilities tool implemented in Terrset and used in the context of CA_Markov allows to consider in the uncertainty of input maps in the quantity of change estimation. However, this introduces important modifications in the calculated quantities of change, resulting in their overestimation ). Accordingly, this method can introduce more uncertainty than the one for which it finds an answer.
Metronamica does not include any method for quantity of changes estimation. The user enters the total number of cells (persistence + changes) that will make up each function class. The tool the user employs to decide the total number of cells will determine the uncertainty of this data. In addition, the transitions modelled for every category will rely  3.3 Allocation of changes LCM, CA_Markov and Metronamica follow a deterministic procedure for change or land use allocation. Therefore, in all executions, the three models allocate the land uses and changes in the same positions (Table 3). LCM and CA_Markov make use of a Multi Objective Land Allocation (MOLA) mechanism. It selects those pixels with the highest potential to change, solving conflicts between different objectives based on the minimumdistance-to-ideal-point-rule (Eastman et al. 1995). factor added to the change potential maps allows transitions to take place in areas less likely to change.
CA_Markov includes a contiguity filter as part of the allocation process, forcing pixels with lower potential values to be simulated as change if located next to previous pixels of the simulated category. In LCM, a zoning layer can be included in the allocation of changes step, multiplying the values of the change potential maps.
The allocation functions of Dinamica EGO (patcher and expander) include a stochastic component to account for the unpredictability of human decision-making. They also include a Cellular Automata component, favouring the simulation of pixels adjacent to land use classes of the same category. If the stochasticity included in the calibration of Dinamica EGO is very high, the user may face problems trying to understand how the model allocates changes and, accordingly, how every parameter influences the obtained simulation. In our exercise, in 20 model executions, only 7.6% of the pixels simulated as change were allocated in the same place, with 34% of the changing pixels allocated in the same place less than 10 times (Table 4 and Fig. 3). The presence of many pixels with similar change potential values may explain this behaviour. Controlling this stochastic component is not easy. It is associated to a specific prune factor, which can be managed by the user, but also to a Monte Carlo approach of land use allocation and the parametrization of the expander and patcher functions (García-Álvarez 2018). To reduce this stochasticity to a minimum, the user must choose a prune factor of 1, define clear and different transition potential values for the candidate cells and parametrize the expander and patcher functions according to the pixel size.
The agreement between simulations across models is not high and depends on the compared models (Table 3). CA_Markov and Metronamica are the models simulating more similar changes. On the contrary, CA_Markov and Dinamica EGO are the models simulating more different changes. These differences cannot be explained by differences between change potential maps, as there is not a direct relation between correlation of change potential maps and the agreement of simulated changes.

Pattern simulation
CA_Markov, Dinamica EGO and Metronamica include a Cellular Automata component to replicate the real LUC pattern. This is lacking in LCM, which however allows to include a dynamic factor of distance to any of the map categories. LCM can infer from this variable the relation between LUC changes and the distance to cells of the other categories. However, this is calculated automatically by the model and therefore dependent on input data uncertainty. Because there is no user intervention possible, there is no direct control of the modelled pattern. In all cases, this CA component or attraction factor, has been the most relevant in explaining the simulated dynamics of Cape Town.
For our modelling exercise, although LCM simulated a general pattern very similar to the reference landscape, the simulation includes many small and scattered patches that do not fit the common pattern associated with land uses like urban residential. Even if we can check a visual coherent pattern in Figure 5, the spatial metrics reveal how LCM was the model that simulated the most new patches of urban residential: 120 opposite to 7 new patches of change between the reference maps of 2002 and 2013 (Table 4).
In CA_Markov the user controls the compactness of the simulated pattern through a userdefined contiguity filter. It up-weights the land use potential values of pixels close to pixels of the considered class and down-weights those which are far from this (Camacho Olmedo and Mas 2018b). It applies the same compactness logic to all modelled classes.
In our study case, urban residential and urban informal changes were simulated according to the same pattern: patches of both classes became more compact, with a reduction in the total number of patches and a bigger mean path size (Table 4).

Metronamica allows to define neighbourhood interactions between all classes of the map
and each function class, making it possible to get a specific pattern for each class. In our model application, urban informal pattern is more scattered (large increment of the number of patches and a lower mean patch size) than the urban residential one. The user can play with the weight of self-attraction rules and the random factor to facilitate the production of new patches in Metronamica. However, this strategy has not always been very effective (Hewitt and Díaz-Pacheco 2017).
In our study case, both CA_Markov and Metronamica faced difficulties when trying to simulate changes as new patches. In the two models, all changes were simulated as infill of existing patches or as an organic halo from them (Fig. 5).
Dinamica EGO simulates the desired pattern through two different functions: expander and patcher, which can be used together or independently. The expander function simulates changes as expansion of previous patches of the same use, whereas the patcher function simulates changes as new patches, disconnected from previous pixels of the same use. The user must indicate to the model the size and shape of the new patches for each modelled transition through three parameters: mean, variance and isometry.
Accordingly, Dinamica EGO is the model that gives more control to the user regarding pattern simulation. For our study, the simulated landscape of Dinamica resembled quite well the reference landscape (Table 4), being maybe the best model when it comes to this point.  Metronamica is the only model that explicitly includes a tool for scenario management.
It also includes a tool to inform about the uncertainty associated to single run simulations.
Dinamica EGO, because of its flexibility, can be designed to produce similar results.
CA_Markov and LCM are more constrained to this end. LCM is only able to simulate business-as-usual scenarios. In CA_Markov, different model applications must be set up to account for different scenarios.

Discussion
Each of the four model software packages compared conceptualized the systems and processes to be modelled in a different way, which resulted in different outputs. These sources of structural uncertainty are discussed in detail in the section 4.1. In addition, all models provide different methods or tools to deal and manage the possible uncertainties of the modelling exercise. They are discussed in section 4.2.

Structural uncertainty
How a system is conceptualized in a model comes with an important source of epistemological uncertainty. CA_Markov, Dinamica EGO and Metronamica approach the phenomena from the same CA based common theory. That is, a complex system's behaviour can be explained by the relation between every of its components (conceptualized as cells) and their neighbourhoods. If this assumption does not lie on the base of the dynamics of the modelled system, these models will probably fail when modelling the intended dynamics. Nevertheless, a wide range of land dynamics have shown to align with this CA behaviour, especially urban sprawl and deforestation (Barredo et al. 2003;White et al. 2015;Kura and Beyene 2020).
In our study case, neighbourhood interactions were the main factor explaining the land use dynamics. However, their importance was different for each category actively modelled. Whereas very important in the case of residential areas, townships and industrial areas, the CA component was less important to explain informal settlements.
Accordingly, models' fit was poorer for this class (Annex 3). When analysing the complexity of CA models regarding the simulated pattern, authors usually focus on the ability of the model to simulate emergent growth isolated from previous developments (García et al. 2011;Hewitt and Díaz-Pacheco 2017). In addition, Pontius Jr. and Malanson (2005)   does not correctly finds this relation, the room for user intervention to avoid this incorrect behaviour is limited. This makes the model very dependent on input uncertainties from data, short calibration periods, etc. For our modelling application, although the modelled result was successful, specific pattern inconsistencies in the obtained simulation could not be corrected because of the model's constrained nature.
Dinamica EGO also incorporates different statistical or automatic methods to find the relation between changes and drivers of change, such as the Weights of Evidence of the Genetic algorithm. However, the user does have the opportunity to change the relation that those methods found to avoid some of those uncertainties. As in our study case, this user intervention may modify to a large extent the model outputs, which will make them very dependent on the uncertainties associated to that user intervention.
In Metronamica and CA_Markov, if following the multicriteria evaluation method recommended by the model's authors, change potential maps mostly rely on user knowledge and understanding, which avoid the uncertainties that come from the method.
However, uncertainties come from user or expert input (Botterweg 1995;Sohl et al. 2016;Li et al. 2017). They may be even bigger than the ones associated to the selected method for change potential calculation.
In automatic or statistical approaches, the model's success will depend a lot on the data and, specifically, on the number and statistical representativity of observed changes. If these are not large or representative enough, the relation that the model finds between changes and drivers of change may be biased. However, this may be a common feature in those models trusting user knowledge if s/he heavily relies on historic dynamics to manually calibrate the models. In these cases, expert judgment could be a solution, as experts can inform about the plausibility of the user's parameters and the simulation results.
According to Botterweg (1995), user's calibrations are only valid for those users or experts who made the calibration, limiting the repeatability of the model application.
However, in our exercise, expert-based calibrations showed high correlation, even with change potential maps obtained through statistical or automatic approaches. On the contrary, change potential maps obtained through different statistical or automatic approaches showed very high variability. These findings are similar to the ones obtained by Krueger et al. (2012). Thus, although the repeatability and easiness of calculation are This manuscript has been submitted for publication in "Environmental Modelling & Software". Subsequent versions of this manuscript, after acceptance and review by the journal, may have slightly different content. If accepted, the final version of this manuscript will be available via the `Peer-reviewed Publication DOI' usually some of the common advantages pointed out to choose automatic or statistical approaches, this is not always the case.
Expert or user knowledge has been pointed out as a useful tool for uncertainty assessment  as far as they accept any map produced through external tools. In Metronamica, although constrained by the specificities of the model, the user also has large flexibility to introduce external maps to calculate the change potential.
LCM supplies three different methods for change potential calculation. However, logistic regression is just provided for pedagogic purposes (Eastman and Toledano 2018b) and developers strongly recommend the Neural Networks option (Eastman 2015). Two extra machine learning techniques have also been provided in the last release of the model (Eastman and Toledano 2018b). However, because machine learning behaves like a black box (Kim 2010;Mozumder et al. 2016), users do not really understand why one method produces different results than the other, which impedes a correct understanding of the model and limits the effective management of the structural uncertainty.
Dinamica EGO, in addition to admitting external maps, also offers different methods to calculate the change potential, which gives the user total control of the structural uncertainty that this step may convey. Although already developed in different studies, CA_Markov and Metronamica do not include alternative options for change potential calculation apart from the MCE and the defined transition formula.
Regarding the quantities of change, the four models allow the user to modify the estimated quantities. He can use this prerogative to employ different methods to calculate those quantities. However, they must adapt in the end to the model requirements: quantities must be provided as a Markov matrix in CA_Markov, Dinamica EGO and LCM and as total number of cells that make up each category in Metronamica. Implementing these methods for change estimation in a dynamic way, connected to the modelling exercise, is very difficult in LCM and Metronamica.
Finally, allocation of changes follows strict procedures in CA_Markov and LCM. The user cannot deal with the structural uncertainty of this step. Metronamica allows the user to change the transition potential formula to deal with this issue, whereas the Dinamica EGO framework is very flexible and can even be used by the user to propose a different allocation of changes algorithm.

Stochasticity as a means to account for the model allocation uncertainty
Stochasticity is considered as an important feature by several authors to replicate real phenomena (García et al. 2011;Renard et al. 2013). It accounts for changes to happen and has been pointed out as useful when communicating the uncertainty of the change allocation step (White et al. 2015). Dinamica EGO is based on a Monte Carlo method that compares the transition probability for each cell with a random number (García et al. 2011), allowing different cells with similar potential to change to be simulated at each model execution.

Process variability uncertainty
At the regional or global levels, the aggregation of uncertain local human decisions brings about new drivers or processes of change, which finally change the foundations of the systems. Accordingly, real systems are far from equilibrium systems, which can evolve to new stages governed by new rules and processes (White et al. 2015). Usually, models deal with this process variability uncertainty by accounting for randomness in the modelling process (Hewitt and Díaz-Pacheco 2017)  Stochasticity facilities those tipping points to happen, allowing to replicate more complex systems. However, changing the foundations of a system is in the reviewed models only possible by means of entering a large randomness, which at the same time introduces great uncertainties in the user comprehension of the model. Neither CA_Markov nor LCM are able to introduce the stochasticity needed to change the system's foundations. The stochastic perturbation approach included in Metronamica also introduces limited stochasticity, as pointed out by Wu (2002) and García et al. (2011) and specifically assessed for this model in our exercise and the study carried out by Hewitt and Díaz-Pacheco (2017). It is important that the model developers provide the users guide and assistance when making this validation exercises. The four software packages include manuals and tutorials that give some tips to this end. However, specific guidelines about validation and, above all, about uncertainty analysis, are lacking. When available, they focus on one or a few analyses and do not make the user aware of the complexity that a full validation and uncertainty analysis may entail.

Communication of uncertainty
There is still a lack of attention in the provision of tools to communicate the uncertainty that the models provide. No model provides enough tools to communicate most of the uncertainties of the analysis to the audience, from the problem conceptualization to the model validation. Models just focus on specific sources of uncertainty, but not for the whole uncertainty of the modelling exercise.
Metronamica includes a tool to obtain probability of change maps based on running the same simulation several times. This allows to account for the stochasticity uncertainty and, therefore, would be especially meaningful in the case of Dinamica EGO, due to the important stochasticity that the model can convey. Notwithstanding, this result could be inconvenient if it is not properly used. It can give the audience a false perspective about the uncertainty of the simulation. It just accounts for the system's uncertainty that the models try to replicate through a stochastic component. However, it does not account for all the other sources of uncertainty which we have addressed in this paper.

Conclusions
We have compared four common LUCC modelling software packages to understand their structural uncertainty and the options that they offer to manage this and other sources of uncertainty. The results showed that each model conceptualized the modelled system in a different way, which led to differences in the way the LUC dynamics and changes were simulated. There is not a best modelling approach, but each model entails different uncertainties and limitations, which must be carefully considered by its potential users.
This manuscript has been submitted for publication in "Environmental Modelling & Software". Subsequent versions of this manuscript, after acceptance and review by the journal, may have slightly different content. If accepted, the final version of this manuscript will be available via the `Peer-reviewed Publication DOI' Statistical or automatic modelling approaches did not prove to be more effective than models relying on user knowledge, neither in terms of model repeatability nor uncertainty management. Models offering more options for change potential calculation, quantity of changes estimation and change allocation procedures were considered preferable, as they give the user more choices to manage the structural uncertainty associated with the system conceptualization that each model encompasses. Of all four software packages, Dinamica EGO stands out regarding these options, followed by Metronamica and CA_Markov.
Randomness and scenario management were identified as two important elements to account for the uncertainty of the modelled processes. However, they are only included in Dinamica EGO and Metronamica. These are also the models offering more options to control the pattern simulation, with Dinamica EGO providing more flexibility to the user to this end.
Despite its constrained nature, LCM captured the dynamics of the modelled system well and delivered a simulation result that resembled the reference map well, while being also very easy to use. CA_Markov did also deliver a result close to the reference. However, this modelling approach was considered too simple regarding different aspects, such as pattern simulation or management of structural uncertainty. Dinamica EGO proved to be an effective tool that, however, if not well managed, may introduce important sources of randomness in the simulation, hampering the understanding of the simulated dynamics.
Metronamica relied a lot on user knowledge, which makes this modelling approach very dependent on the quality of the user calibration.
All simulation results were different according to the chosen modelling approach, both in the allocation of simulated changes as in the production of change potential maps. These differences may be considered a way to better understand the uncertainties that are usually associated to a specific LUCC modelling exercise, given the structural uncertainty of each specific LUCC modelling software.
In all cases, we have identified a lack of attention to important aspects related with uncertainty management, such as the communication of model uncertainties and the provision of tools and guidance for uncertainty analysis.
This analysis is limited by the specific study case selected for the model comparison and its specificalities. In this regard, model outputs were only compared for one specific case and one historic period. Further research should look at the differences between model This manuscript has been submitted for publication in "Environmental Modelling & Software". Subsequent versions of this manuscript, after acceptance and review by the journal, may have slightly different content. If accepted, the final version of this manuscript will be available via the `Peer-reviewed Publication DOI' executions when using different historic periods and study areas. Assessing the plausibility of model parameters and results, based on expert judgment or other strategies, could be also an interesting point to assess the uncertainty that each model software package entails, as we have only judged the models' success based on their quantitative performance with respect to a historic period of reference.