In the past, the first concern about agricultural land was its potential crop productivity and this issue was dealt with land evaluation. In 1976 the Food and Agriculture Organisation (FAO) defined land evaluation as (FAO, 1976): the assessment of land performance when used for specified purposes. It involves the execution and interpretation of basic surveys of climate, soils, vegetation and other aspects of land in terms of the requirements of alternative forms of land use.
Since then, land evaluation has taken on a leading role in the sustainable management of natural resources. In fact, it integrates information from different fields and at different spatial scales, including: type of soil, climate data, vegetation, geomorphology, economic and social drivers, to find an optimal balance between the needs of human societies and the preservation of natural resources in the long term. From land planning to evaluate specific uses of the land, to the improvement of their management and the prevention of soil degradation dynamics, land evaluation over the years has become a useful tool to ensure a sustainable exploitation of a territory and to preserve its productive and environmental functions (Giordano, 1999; Eliasson, 2007). Specifically, a rational land evaluation should provide an initial analysis of the biophysical component of the territory, to determine the suitability of the soil in relation to generic uses, such as agricultural, forest or natural uses (land capability) or to specific uses (land suitability), followed by a second step defining possible alternative uses in consideration of different socio-economic scenarios (Giordano, 1999).
Land evaluation deals with issues to which there is no unique answer. The dynamism of human societies and market conditions, together with the variability of environmental features, at different scales, set limits to our ability to predict specific use destination for a given territory. Moreover, our spatial knowledge of the territory is not perfect and the temporal dimension further expands the margins of error in any evaluation. In addition, the concepts of use and suitability cannot be quantified in absolute terms, while they are usually referred to specific purposes and developed in comparative terms (Giordano, 1999). Therefore, over time, the increasing complexity of such an analysis has set the limits of conventional approaches, while scientists have started to investigate climate change and land use and cover change phenomena. From the initial focus local productivity and rational land use concerns, scientists’ efforts have turned towards the study of complex systems and global dynamics, thus developing new methods and tools. Among these, the agro-ecological zoning (AEZ) methodology was developed by FAO in collaboration with the International Institute for Applied Systems Analysis (IIASA). The project dates back to 1978, but only in the last twenty years this assessment methodology has spread around the world, thanks to the potential offered by the latest technologies and, in particular, by spatial databases of resources available on a global scale (Eliasson, 2007). The AEZ methodology consists of separating a territory into spatial units of small size and homogeneous within them for soil suitability compared to specific uses, environmental impact and potential production (Eliasson, 2007).
The recent technological improvements have also allowed the development of software tools to monitor the condition of current resources and the evolution of land use systems (Verburg and Van De Steeg, 2009). In this context, modelling has become a tool to explore different scenarios and alternatives, to improve decision-making processes in land use, land cover changes and global change management (Verburg et al., 2006). Therefore, the use of mathematical models have allowed better understanding of the drivers that determine these events and scenario analysis has become a widespread tool for exploring land use systems and their potential evolutions. Scenarios are not used to make predictions but rather to direct the decision making process in agreement with different alternatives, in contexts where uncertainties are managed according to specific assumptions. Scenarios then allow the comparison of potential consequences that would result from decisions regarding land use, land cover changes and global change management (Verburg et al., 2006).
Several models have been developed. Most of them use low-resolution data (>1×1 km), because of the difficulty in obtaining high resolution ones. These are models useful to build land use scenarios from global to national scale. Few are instead the models that employ data from high-resolution raster maps, with a pixel size lower than 1x1 km; among these, probably the best known is CLUE-S. Unlike CLUE, useful from global to continental scale, the CLUE-S is not based on sample surveys and census data but on raster maps stored in a geographic information system (GIS) environment with a hierarchical structure, to analyse the evolution of land use systems with a multi-scale approach (Verburg et al., 2002). The CLUE-S model calculates the probability for each pixel to belong to a defined land use class, according to a logistic regression in which the independent demographic, economic, technological, political, institutional, cultural and biophysical variables are the drivers of change (Verburg et al., 2002). For instance, the allocation of paddies in land use systems is modelled according to elevation, slope, distance from cities and demographic density, but without consideration of soil properties.
Another model that employs databases with a resolution lower than 1×1 km is the Patuxent landscape model. This model focuses on the reconstruction of the ecosystemic component (ecological module) in relation to economic and demographic factors (economic module), to guarantee the best land use management to preserve good water quality, at basin scale. Within the ecological module, the model includes information such as ecological successions, quantity and quality of surface water, the depth of aquifers, health conditions and vegetation growth, health conditions of habitats and their fragmentation. Economic and demographic dynamics concerning land uses and agricultural practices are instead processed within the economic module (Voinov et al., 1999). The geomorphologic and biophysical informational component of soils is not included in the model.
At global scale, some modelling efforts have been recently published, but the use of a low spatial resolution database (from ~10 × 10 km to ~55 × 55 km) affects results because of excessive approximation. In fact, these models are not oriented to the modelling of land use systems in detail, but only of the main environmental dynamics at global scale (Meiyappan et al., 2014).
Taking all this into account, both AEZ methodology and model-based approaches appear as useful tools to analyse the impacts of climate change, land use and cover change phenomena on socio-ecosystems and on global biochemistry at different spatial and temporal scale. However, it appears also that most of recent studies do not consider soil features at large scale as well as biophysical and geomorphological features, thus loosing the possibility of considering some of the most important phenomena, which determine the capability of biological systems to provide mankind with resources. This is often due to the lack of high-resolution data, but more broadly it is due to the limited efforts put on the integration of multidisciplinary knowledge. A likely consequence is a bias in the results obtained, affected by the lack of consideration of important factors - erosion phenomena, for instance - limiting potential production, in particular when land use dynamics are analysed over the long time periods typical of climate change studies. For example, most of socio-economic studied on agriculture and climate change studies do not consider any long-term change in land suitability and productivity and their feedback effects on socio-ecosystems. On the contrary, the European Thematic Strategy for Soil Protection whose related proposal for a common framework directive [COM(2006)232], identifies erosion, soil organic matter decline, compaction, salinisation and landslide events as the main risks to soil conservation and long term productivity of agro ecosystems (European Commission, 2006).
With the aim of overcoming some of the problems mentioned above, in this work we explored the possibility of developing virtual territories with suitable software, statistically replicating the morphological, ecological and land cover features of real environments, on which to explore scenarios of land use change and climate change phenomena. For greater clarity, we intended the statistical replication of real environment features as a way to generate similar features in a new, not existing and virtually developed system, the so-called virtual territory. We define it as an infrastructure of coherent data with internal consistency between spatial variables, to which we can apply existing models or test new ones in order to carry out scenario analyses. The ambition is to use such models in a controlled spatial environment as an analogue of punctual agronomic experiments in parcels or laboratory. As in the case of punctual experiments, there is no will to substitute studies on real world territories, but just to complement them with more opportunities to explore the complexity of agro-ecosystems and provide opportunities for what-if simulations. In order to guarantee a robust statistical link between the real and the virtual territories, the high-resolution databases of biophysical and geomorphological information available in well-studied areas should be studied to assess the relationships among the most important variables.
Provided that strong statistical relationships between spatial variables can be described, they can be later on used to guide the design of internally consistent virtual analogues and of variants of existing landscapes, which can be utilised for scenario analysis. Therefore, the proposed approach is intended to expand the opportunities for scenario analysis, overcoming the constrains mentioned above limiting the possibilities to conduct long-term analyses in the real world in fields related to global change sciences, such as climate change adaptation in the agricultural sector.
In the Materials and methods section we first present fractal theory as the method to build a virtual digital elevation model (DEM) as the geomorphological basis of the virtual territory and we describe the area selected as a reference for the application. Later on, we introduce the methods adopted to build the GIS internally consistent layers with all the information required to generate a virtual dataset with statistical fitting to those of the existing territory and we introduce the erosion model adopted to demonstrate the approach. In the Results section we present the scenarios considered and the outputs of spatial simulations as raster maps previously generated. The proposed approach and the results obtained are discussed in the fourth section of the paper, followed by some conclusive remarks and the identification of future research needs.
Materials and methods
Fractals and the real world
Real world objects, such as a landscape identified as a portion of land surface, may exhibit two different complexity types: fractal and not fractal. Non-fractal complexity is the result of a series of separated and uncorrelated events that have occurred over time. Because of the mutual independence of these events, it is really difficult to interpret the features that have determined this complexity. Conversely, fractal complexity characterises those features that exhibit statistical self-similarity, namely similar shapes repeated within a specific range of scales, defined as complex geometric objects resulting from the repetition of the same morphology at different scale of magnitude, between a lower and an upper crossover scale (Ebert et al., 2003). A wide variety of natural phenomena show fractal structures, such as for example: lightning, trees branches, and vascular systems, but also the morphology of natural surfaces (e.g., mountains) are known as statistical self-similar fractals (Ebert et al., 2003) and their local fractal dimension can be estimated using statistical calculation procedures (Eastman, 2009). This is estimation because a real surface is the result of a combination of different events, which, as stated above, may determine both fractal and non-fractal shapes. If a natural landscape can be described through the parameters of its fractal dimensions, then the reproduction of those shapes and phenomena should be possible through the procedural generation of synthetic landscapes, based on algorithms that generate fractal surfaces. In practice, the shapes of a fractal surface depend on the basis function, a well-defined mathematical function identified by a specific amplitude and spatial frequency (Ebert et al., 2003). The fractal dimension, also called roughness, is the main property that characterises a fractal surface and it is given by the Euclidean dimension 2 to which a fractional part is added (e.g., 0.3), defined as fractal increment. Procedural fractal surface generation consists of an iterative process where the shape defined by a basis function is repeatedly added and rescaled to build a complex surface. The amount with which spatial frequency changes between subsequent scales are defined by lacunarity: with a lacunarity set at 2, the spatial frequency doubles while the added shape size halves (Ebert et al., 2003).
Real topographic surfaces exhibit a variable local fractal dimension because of their heterogeneity deriving from geological features and other morphogenetic phenomena (Gallant et al., 1994).
Therefore, multiple fractal dimensions are required to generate virtual heterogeneous surfaces as in the real world. Musgrave pioneered the development of procedures and computer codes to build first multifractal models in order to increase realism, with the aim of reproducing erosion dynamics with so-called ridge multifractal models (Ebert et al., 2003). The present work develops upon those seminal works, and adopts a multi-fractal approach to generate realistic virtual territories.
Methods and tools for the generation of fractal surfaces to reproduce morphological features of an existing territory
Having opted for a methodological approach based upon multifractals we reviewed the existing software tools, but we did not find a solution allowing us to have full transparency and control of the implemented procedure and providing all the elaboration routines required. To build the virtual morphological surface, i.e., DEM as defined in GIS technology we adopted the Blender software and its plug-in ANT Landscape, while the other required spatial data processing functions were found in two GIS software tools: QGIS with its GRASS extension, and TerrSet.
Firstly, a portion of territory of the Veneto Region was identified and its morphological features quantitatively analysed. Secondly, a series of virtual DEMs ware generated by using the ANT Landscape tool with varying parameters, and then compared ex post with the real world, focusing in particular on their fractal dimension values and related indices. The virtual DEMs, along with the DEM of the existing area, were decomposed in the frequency domain with the 2-dimension-fast Fourier transform (2D-FFT) routine provided by the TerrSet GIS software. The 2D-FFT is a useful tool for analysing the hierarchical structure of territory components of different scale (Florinsky, 2012), by decomposing the cosine (real) and the sine (imaginary) components. The module returns a cosine component raster image, a sine component raster image and a power spectrum raster image (Eastman, 2009). The power spectrum image can be interpreted in the following way: i) pixel values represent power spectrum according to the following equation:2009); consequently, assuming that these features are shared between two DEMs, their power spectrum identifies the same frequencies.
The power spectrum is interpreted as the degree of resonance of all different frequencies (Eastman, 2009), namely the hierarchical structure of the geomorphological components of the DEM, no matter what the original position of each shape is. Low frequencies identify smallscale morphologies and average slopes while high frequencies are interpreted as high scale shapes and higher slope values (Toth et al., 2014).
Having identified the virtual DEM with the best fit to the original landscape, we moved to demonstrate its use by implementing a simple model for erosion simulation over long time period and under different climatic scenarios. AEZ were adopted to provide a concise description of the combinations between the most important physical features of both the real and the virtual territories and allow for cross-comparisons.
Development of scenario simulations on the virtual territory
The most widely used empirical model to quantitatively estimate soil erosion is the universal soil loss equation (USLE) (Wischmeier and Smith, 1978), and its more recent evolution into the revised USLE (RUSLE) (Renard et al., 1997). Both those models estimate the annual loss of soil, caused by sheet and rill erosion but the RUSLE was specifically designed to exploit the potential of GIS software for spatial analysis (Bosco et al., 2008) and it was thus selected for this study:
A = annual soil loss (t▪ha–1▪yr–1)
R = rainfall erosivity factor (Mj▪mm▪ha–1▪h–1 ▪ yr–1)
K = soil erodibility factor (t▪ha▪h▪ha–1▪Mj–1 ▪ mm–1)
LS = slope length factor (dimensionless);
C = cover management factor (dimensionless);
P = human practices aimed at erosion control (dimensionless).
Effects of hypothetical climate changes scenarios were implemented by recalculating the R factor, assuming monthly variations of available rainfall data for the existing territory to reproduce the effects of an exacerbation of the characteristics of Mediterranean climate, with intense precipitation concentrated in shorted periods over the autumn and winter seasons, as foreseen in the most recent report of the Intergovernmental Panel for Climate Change (IPCC, 2013).
Data stored in the digital soil map of the Veneto Region, with detailed information from regional to sub-system level were used to describe land suitability through the identification of AEZ for rainfed agriculture, following the criteria proposed by Eliasson et al. (2007). AEZs were also adopted as a key to provide a concise description of land features and to maintain internal consistency between the variables of the virtual territory. Moreover, AEZs were used to monitor changes in land suitability over time, as affected by the simulated erosion phenomena. Results were stored in a vector GIS layer with polygons characterised by unique combinations of six AEZ criteria (slope, soil depth, soil fertility, soil drainage, soil texture and soil chemistry), identified by a 6-digit code (e.g., 111111 for the AEZ with optimal characteristics, or no limitations).
Selection of an existing territory and generation of statistically similar virtual territories
An existing territory has to be identified and studied in order to provide the information needed to generate a coherent virtual dataset with physical features statistically similar to the real one.
We selected an existing territory, located between the landscape areas of Monti della Lessinia and Alta Pianura Veronese included in the Piano Territoriale di Coordinamento Regionale (PTRC) of the Veneto Region, with the aim of using regional to local datasets. Figure 1 shows the selected area (Regione Veneto, 2013).
The area exhibits a positive altitude gradient from south to north, until reaching about 1549 m. Alpine foothills (Prealpi) slope down towards the south plains around Verona with significant changes in vegetation and land uses in relation to elevation. There is a high density of small urban cores in the valleys, under 1000 m of altitude. In the plain rural areas are opposed to more densely populated urban centres around the city of Verona. The geomorphology is mainly affected by tectonic movements and the erosion of rainwater. The selected area is generally constituted by calcareous lithotypes, which determine the almost total absence of surface hydrography (Regione Veneto, 2013).
By using the parameters in Table 1, we generated six virtual DEMs with varying values for the parameter Dimension.
Linear regression between the fractal dimension and the average values of the power spectrum calculated on raster maps of virtual DEMs with 2D-FFT was calculated to guide the generation of the multi-fractal DEM (Figure 2).
According to this statistical relationship, we calculated an average power spectrum of 15.04 on the DEM of the selected area and we estimated a Dimension of 0.41. With the highest fractal dimension of 2.59, we generated a new virtual DEM in ANT landscape. After having georeferenced it and rescaled the vertical scale to the elevations of the selected area, we got a virtual DEM, with the statistics comparable to the real one (Table 2 and Figure 3).
With the module r.regression.line of the GRASS plug-in of QGIS, we estimated a Pearson’s correlation coefficient of 0.79 between the power spectrum maps of the existing and of the virtual DEMs.
Implementation of the revised universal soil loss equation model and agro-ecological zoning classification
R factor defines the influence of climate on erosion phenomena, combining the effect of the mechanical action of rainfall with superficial runoff, both laminar and rill (Bosco et al., 2008). As previously done by Borin and Bonamano (2005) to estimate the rainfall erosivity factor for some areas of the Euganean Hills (Veneto Region, not far from the selected area), the R factor Mj▪mm▪ha–1▪h–1 ▪ yr–1 was calculated by using the following method proposed by Arnoldus (1980) on monthly rainfall data, available for the selected area from 1994 to 2012:
This equation was useful to describe rainfall seasonality patterns both in the current climate and in the hypothetical ones generated for scenario. We estimated R-values of 287.10 Mj▪mm▪ha–1▪h–1 ▪ yr–1 for plains and hilly areas and 668.22 for mountainous areas of the selected territory. Those values were then applied to corresponding areas of the virtual landscapes with the raster calculator of QGIS, thus producing the R factor map of the virtual DEM.
LS factor represents the area where the superficial runoff diverges or converges and defines the upslope contributing area estimated according to Bosco et al. (2008) with the following equation:
Flow_acc = flow accumulation (the number of pixels through which water flows towards a specific pixel);
a = slope;
cell_size = pixel size.
Starting from the virtual DEM, the flow accumulation map was calculated by using GRASS module r.watershed, while the slope map was generated with the geomorphological analysis tool of QGIS, and the LS factor map was generated by using the raster calculator.
C factor is the ratio between soil loss in defined land cover conditions and the erosion that would occur on a soil without any protection. To generate the C factor virtual map it was necessary to reproduce a coherent land cover distribution starting from the CORINE land cover (CLC) 2006 map for the selected area, with the aim of assigning C factor values available from literature (Bosco et al., 2008) to land cover classes. Having produced the required set of map layers as described above, the map of AEZ of the virtual territory was produced by means of map overlaying, with the coding system presented in Table 3, ending up with 33 different AEZs. Afterwards we performed a cross-tabulation between the AEZ raster map and the CLC map reclassified in two levels, with the aim of calculating the percentage proportion of CLC category distribution. The cross-tabulation of AEZ versus land cover classes for the existing territory was used as a key to reproduce realistic combinations between those two layers in the virtual landscapes, in order to assign meaningful C factor values. Therefore, we created a random distribution of 1000 vector points within the extent of the virtual DEM and with such points we generated a layer of Voronoi polygons. Afterwards, we overlaid that layer on slope classes used for AEZ classification and we assigned the slope class value to each polygon of the virtual DEM and we used the resulting map to randomly distribute AEZ codes, with the constrain of meeting their proportions in the real world. In this way we obtained a vector layer that shows a new virtual distribution of AEZs, consistent with the slopes of the virtual DEM. Eventually, the cross-tabulation between AEZ and land cover guided the assignment of CLC categories to the virtual territory and the calculation C factors.
The K factor was calculated on the basis of soil properties derived from the digital soil map and assigned through the virtual map of AEZ. It was estimated according to Romkens et al. (1986) on the basis of a regression conducted on a large global dataset of K factor values, which produced the following equation [revised by Renard et al. (1997)]:
Dg is the weighted medium diameter compared to dominant texture classes (mm);
di is the highest diameter of each texture class (sand, silt and clay) (mm);
di-i is the lowest diameter of each texture class (mm);
fi is fraction of each texture class. With the equations (5) and (6), we estimated a K factor value of 0.0088 for AEZs which have soils with coarse texture and 0.0404 for AEZs with medium to fine texture.
In this study we did not entered into the detail of considering different management practices and then the P factor was set to 1 for all AEZs. The raster calculator tool of QGIS was used to run the calculation of erosion, producing a RUSLE map, which presents the average loss of soil in tons per hectare per year for the virtual territory (Figure 4). The results were compared to those obtained by ARPAV (2008; P.29) for the reference territory, obtaining a good match, with values for the plain below 2, and increasing values in steep slope area, which could episodically fall into the class above 40 tons.
The final erosion calculation concerned the estimation of the depth of eroded soil, in order to assess whether erosion phenomena could determine a change in the AEZ classification. Assuming a bulk density of 1.38 g▪cm–3 for the AEZs with medium to fine texture and 1.61 g▪cm–3 for those with coarser textures we calculated the tons per centimetre of depth on each pixel area: 11.62 for the AEZs with finer texture; and 13.55 for the AEZs with coarser texture. Figure 5 presents the spatial results of the simulations. The first map presents the initial soil depths of the virtual territory, while the maps below show the results of the RUSLE model over a time span of 50 and 100 years with the current climate. Climatic scenarios were analysed by varying monthly rainfall data as mentioned above, to apply Arnoldus’s method (1980) to recalculate R values representing dryer summers and wetter winters, together with limited decreases of annual rainfall (Table 4).
Table 5 shows the relative and the absolute amount of soil eroded after 100 years of erosive processes, in the different rainfall scenarios.
While simulations with current climate showed only very limited effects on AEZ classification caused by variations of soil depths as a consequence of rainfall after 100 years, hypothesis 1 and even more hypothesis 2 climates showed non negligible effects. Current climate variations of AEZ surface areas ranged between 0 and 9.3%; hypothesis 1 generated changes in soil depth causing changes in AEZ areas up to 12.2%; and hypothesis 2 produced changes up to 24.3%. Running the RUSLE model on the virtual territory for 100 years, it turns out that 14.01% of the total extent is still assigned to the most suitable AEZ for rainfed agriculture code 111111. On the other hand, 0.23% of the surface is completely eroded, while many AEZs undergo a transition in terms of soil depth categories. AEZ 211111, which has only a slight slope constraint, undergoes the highest reduction of 9.3%. Other significant reductions are 3.7% for AEZ 112122, 3.3% for AEZ 321111 and 3.0% for AEZ 321211. Simulations for hypothesis 1 rainfall showed an increase of 30.44% of the extent of total eroded soils and of 9.02% of the overall amount of tons of eroded soils. In the second hypothesis of rainfall change these increments are respectively 182.61% and 54.50%. Compared to the first simulation, the two scenarios demonstrate how climate evolution hypothesised by the IPCC could increase erosion phenomena, even with reduced annual precipitation.
Discussion and conclusions
Global change studies are research fields that involve multi-disciplinary knowledge and require continuous development of new tools and technologies, to enable us to improve our ability to monitor these phenomena throughout time and explore them in future projections. In this context, we explored the possibility of developing virtual territories, with the aim of statistically replicating the morphological, ecological and land cover features of real environments, on which to test and experiment land use change and climate change phenomena.
The first step in building a virtual territory is the generation of a geomorphological surface statistically similar to a selected existing territory whose features the scientist wants to reproduce. We adopted the solution of procedural fractal surface generation because natural geomorphological shapes exhibit fractal morphologies with statistical selfsimilarity. The limit of this method is that it is not able to reproduce non-fractal morphologies resulting from a series of separated and uncorrelated events.
The virtual DEM we obtained with the highest fractal dimension of 2.59 shows a Pearson’s correlation coefficient of 0.79 between the power spectrum maps of virtual and existing DEMs, which we considered as a satisfactory correlation. At the moment, the procedure consists in a manual repetition of trial-and-error approach to generate virtual DEMs on the basis of fractal parameters calculated on the existing territory, followed by statistical tests of goodness of fit. A more efficient approach would consist in the development of an optimisation procedure based upon automatic sequential generation of DEMs and statistical tests, with defined ranges of fractal parameters.
The approach described above used the RUSLE model as an example to demonstrate how to develop virtual landscapes meeting the spatial variability of the various physical variables and with internal consistency between them, i.e., the between the values stored in the package of GIS layers (soil characteristics, slope, land cover, etc.).
The adoption of AEZs as a means to refer the characteristics of the territories to land suitability and, in case, later on also to its productivity, appeared as a crucial element to provide a meaningful and concise representation of the landscapes, to maintain and assess their realism and to monitor possible changes of agro-ecosystem features over long time frames.
For this study, GIS software showed to be an indispensable tool for the generation and the development of the informational component of virtual territories, even though we could not find all the required procedures within the same package. The development of automated multifractal DEM generation tools as GIS plug-in appears as a stimulating perspective for future research efforts in this field.
Climate change studies are affected by a variety of uncertainty factors. Many attempts have been performed to simulate the behaviour of global socio-ecosystems under global change, but quite often, the ambition to consider all the exogenous drivers and the endogenous dynamics ends up with a substantial loss of control of simulation conditions. An approach based upon virtual territories has opposite pros and cons: it loses indeed the immediate reference to existing territories, but at the same time it allows a much stronger control of the territorial variables and thus it facilitates the exploration of causal chains in systems’ behaviour and thus the analysis of what-it scenarios in controlled environments.
Indeed, the package of GIS layers which constitute a virtual landscape are the result of a series of approximations applied with the aim of reproducing dominant features of an existing territory, but we believe that the proposed approach showed a good potential to complement studies conducted on existing territories. We may in particular imagine the identification of areas of greatest interest at global level, upon which one could develop virtual spatially explicit laboratories to run long-term controlled simulation experiments with varying combinations of environmental (e.g., climatic), socio-economic, and policy scenarios, and their results could significantly consolidate and improve global scale simulations run by existing integrated assessment models.