Modeling the spatial distribution of crop cultivated areas at a large regional scale combining system dynamics and a modified Dyna-CLUE: A case from Iran

Iman Mesgari

Iran University of Science and Technology, Dept. of Industrial Engineering. Narmak-Tehran, 1684613114 Iran.

Mohammad Saeed Jabalameli

Iran University of Science and Technology, Dept. of Industrial Engineering. Narmak-Tehran, 1684613114 Iran.



Agricultural land use pattern is affected by many factors at different scales and effects that are separated by time and space. This will lead to simulation models that optimize or project the cropping pattern changes and incorporate complexities in terms of details and dynamics. Combining System Dynamics (SD) and a modified Conversion of Land Use and its Effects (CLUE) modelling framework, this paper suggests a new dynamic approach for assessing the demand of different crops at country-level and for predicting the spatial distribution of cultivated areas at provincial scale. As example, a case study is presented for Iran, where we have simulated a scenario of future cropping pattern changes during 2015-2040. The results indicated a change in the spatial distribution of cultivated areas during the next years. An increase in the proportion of rice is expected in northern Iran, whereas the proportion of wheat is increasing in the mountainous western areas. Wheat and barley crops are expected to become dominant within the cropping system throughout the country regions.

Additional keywords: integrated modeling; land use; cropping pattern; system dynamics.

Abbreviations used: CA (competitive advantage); CLA (cultivated land area); GDP (gross domestic production); RCE (relative change elasticity); RSA (relative suitability advantageous); SD (system dynamics); TRC (total relative capability).

Authors' contributions: Conceived and designed the experiments: IM and MSJ. Performed the experiments, analysed the data and wrote the paper: IM. Critical revision of the manuscript: MSJ.

Citation: Mesgari, I.; Jabalameli, M. S. (2017). Modeling the spatial distribution of crop cultivated areas at a large regional scale combining system dynamics and a modified Dyna-CLUE: A case from Iran. Spanish Journal of Agricultural Research, Volume 15, Issue 4, e0211.

Received: 23 Oct 2016. Accepted: 28 Dec 2017.

Copyright © 2017 INIA. This is an open access article distributed under the terms of the Creative Commons Attribution (CC-by) Spain 3.0 License.

Funding: The author(s) received no specific funding for this work.

Competing interests: The authors have declared that no competing interests exist

Correspondence should be addressed to Oscar Iman Mesgari:





Material and methods





Agricultural development is the basis of human survival and thus, an appropriate insight into the future changes in the cultivated areas is significantly important for both whole society and policy makers. More specific in agricultural systems, an appropriate insight into the future spatial and temporal arrangements of cultivated areas have been highlighted as a crucial parameter, reaching to which requires scientific methods and models (Castellazzi et al., 2010).

The simulation of land use changes has become a widely used technique in studies of land use planning and ex-ante assessment of new technologies, policy interventions and climate change on agricultural systems (e.g. Doglioni et al., 2009; Zhong et al., 2009; Wang et al., 2010). Some of these models, such as IMPACT (Rosegrant et al., 2008), focus on economic aspects and try to represent a consistent framework for market and trade mechanism, while some others, such as CLUE (Verburg et al., 2002), are spatially explicit and generate the land use patterns based on geographical characteristics. The CLUE model is now one of the most widely applied models with approximately 30 applications spread over the different regions of the globe, addressing a wide range of land-use change trajectories including agricultural intensification, deforestation, land abandonment and urbanization. This model is a tool to better understand the processes that determine the changes in the spatial pattern of land use and to explore possible future changes in land use (Verburg & Overmars, 2007). Beside these improvements, a consensus has currently grown among model developers who identifying that “optimal” or “probable” crop lands arrangements are usually an interdisciplinary activity. This requires to ideally give a balanced representation of the economic, geographical, as well as policy factors over a range of spatial and temporal scales (Parker et al., 2002; Schaldach & Alcamo, 2006; Schreinemachers & Berger, 2011).

In spite of advantages of integrated models, they incorporate the complexities. Senge (1994) discusses two types of complexity: (1) detail and (2) dynamic. Detail complexity is associated to systems with many components. Dynamic complexity, however, is associated to systems with effects that are separated by time and/or space. Complexity of details is associated to agricultural systems because they have many components at different scales and require high data preparation and processing, especially for fine resolutions. On the other hand, complexity of the dynamics is related to this system; having effects which are separated by time and/or space (Zarghami & Akbariyeh, 2012). Due to these complexities, use of many crop studies is faced with limitations for national policy makers who need to understand the consequences of policies and to rapidly test a wide range of scenarios. Dynamic modeling, especially system dynamics (SD), is an effective approach to characterize the whole system by integrating various processes and reflecting the interactions between structures and behaviors of complex systems (Zhang et al., 2002; Mesgari et al., 2017). SD method, proposed by Professor J. Forrester (1958) at MIT, is a powerful tool for modelling and analysis of complex systems that are composed of interacting subsystems or factors that work together to influence overall system behavior (Zhou et al., 2016). SD has been widely used to evaluate different policies/strategies in order to improve system performance (Sharon et al., 2011). The main advantage of SD approach is to describe this kind of systems with high level of uncertainty, causal ambiguity, and structural complexity (Peng et al., 2014).

There are some dynamic simulation models on cropping pattern changes in recent years that tried to combine economic and environmental components at fine resolution (e.g. Luo et al, 2010; Britz et al., 2011; Akıncı et al., 2013; Pilehforooshha et al., 2014). But there is an emerging need to introduce new models at large country-level or provincial scale (Priya & Shibasaki, 2001).

Therefore, this study mainly contributes by describing a dynamic and integrated approach. This approach aims to allow the rapid generation and testing a range of possible scenarios at country-level to provincial scale. Generally, the suggested model can be differentiated into two scales, national and provincial. At the national scale, SD model and national agricultural statistics that are aggregated at administrative levels are used to estimate the cultivated land demand for each crop. The SD technique enables us to design and test different scenarios based on ecological and socio-economic deriving forces, including population growth, economic development, climatic changes, policy rules, and so on. At the provincial or regional scale, a modified version of Dyna-CLUE (Verburg & Overmars, 2009) is designed to simulate the spatial arrangement of estimated land demand between provinces or states of the country.

This paper describes the modelling approach, the inputs, outputs and processes used within the suggested approach. The model evaluation is presented through the analysis of its fitness for this purpose. Application of the model is then demonstrated with a case study on the Iranian agricultural system.

Material and methodsTop

Data and study area

Iran is the second largest country in the west of Asia, with a population of approx. 80 million people. This county comprises 31 provinces and has a total area of 1.648 million square kilometers, of which 8.7% is cultivated land. Average annual precipitation is 228 mm and approx. 90% of the country is arid or semi-arid. In 2014, main crops included wheat, barley, rice and maize, which were grown on 53%, 16%, 5% and 2% of the country’s cultivated areas, respectively. For empirical study in this work, Iran’s 31 provinces were considered. In each province, cropping pattern of four dominant crops (wheat, barley, rice, and maize) were simulated and analysed for the next 25 years. Most of the farms in Iran are monoculture and in-rotation cultivated areas are negligible, then in this study we assumed that all agricultural areas will be cultivated monoculture.

A large amount of data, including national and provincial data, was used in this research (Table 1). Three types of data have been used: (1) some data were input directly into the SD sub-model for demand estimation; (2) some data were used by the crop allocation sub-model to simulate the distribution of estimated demand within the provinces; and (3) some data were used only for calibration of the model.

Table 1. Datasets used in this study

Structure of integrated model

The integrated model consists of two sub-models: SD sub-model and crop allocation sub-model. SD estimates the demanded area for agricultural land at national level by considering the interaction between agronomics driving forces and without spatial resolution. Crop allocation sub-model simulates distribution of demand within geographical regions of the country by considering spatial characteristics of each region. Fig.1 shows the overall structure of the model. The exogenous variables of the model include population growth, economic growth, climate change, and ecological features of each region. The output of the model is presented in two steps. First at country scale, the model estimates the demand, the production, the trade and the area of farmland of each crop. Second at regional scale, the model determines the share of each crop from the total cultivated area of each region.

Figure 1. General structure of integrated model. GDP: gross domestic production. SD: system dynamics.

SD sub-model

Sterman (2000) developed five steps to create SD models: 1) Problem articulation; 2) Dynamic hypothesis; 3) Formulation and deployment in software; 4) Testing; 5) Policy formulation and evaluation.

Recently, Mesgari et al. (2017) have presented a SD model to predict the total agricultural land use demand at the national scale. In this study, a modified version of that model has been used to estimate the demand, production, trade, and cultivated area of each crop at national scale. Fig. 2 describes the causal diagram of this SD model. The major assumptions in the development of this sub-model include: (1) demand for agricultural products will increase as population growth and income growth; (2) increasing demand will stimulate the supply side through the price mechanism; (3) part of the demand will be satisfied through the domestic production and the rest is met by trade mechanisms; (4) to achieve equilibrium, policy makers may apply optimization policies to interact between demand and supply such as: import tariffs, export incentives, and subsidies.

Figure 2. Logical framework of system dynamics model (color variables are exogenous). +: incremental effect, -: reducing effect.

The stock-flow diagram (SFD) is the core of SD model, and is the process of quantization and materialization of the causal loop diagram (Li et al., 2012). Fig. 3 shows the SFD of the model which has been depicted in the Vensim software. Three main feedback loops could be distinguished in this diagram:

(1) Demand loop: In this loop, increasing “GDP” and increasing “population” lead to higher “demand” for crops. GDP and population are exogenous variables and are the main drivers of demand. Consequently, increasing demand cause an increase in the “price” of crops. Then the increased price will reduce demand and balance it in a negative loop.

(2) Supply loop: In this loop, rising price stimulate and increase “production” in a positive loop due to rising “value added per worker” and enhancing “cultivated land area” consequently. Increasing production will reduce price and smooth it.

(3) Trade loop: In this loop, a drop in price (due to decreased demand or increased production) will reduce “ratio of in/out price” and then, import will be economically justified. Increasing import will reduce price and smooth it. On the contrary, this mechanism would be applied to export. In fact, trade in agricultural products creates equilibrium between prices within the country and abroad.

Figure 3. Vensim stock-flow diagram of system dynamics model

The main output of SD sub-model is the demanded area for each crop at country level. This output is the input of crop allocation sub-model. Complete details about the variables and formulas are available in Table S1 [suppl].

Crop allocation sub-model

By inspiration from Dyna-CLUE (Verburg & Overmars, 2009), an algorithm has been suggested to project the share of each region from national demanded cultivated area for each crop (referred in this paper as “crop allocation”). Dyna-CLUE allocates a land use to each cell based on competition between different land uses while crop allocation procedure determines portion of each region area which will be dedicate to a specific crop cultivation based on competition between different regions. In this algorithm, as depicted in Fig. 4, first we must calculate the total available land for each crop in each region based on total agriculturally usable area, spatial policies, allowable conversions and other specific restrictions. Then we must calculate the relative capability or advantage of each region for the cultivation of each crop. Finally, algorithm distributes the demanded area that is derived from SD sub-model among regions taking into account the total available land and the relative capabilities of each region.

Figure 4. Algorithm of crop allocation sub-model

More specific, the crop allocation sub-model like any other allocation mechanism needs to consider the following three parts: (1) cultivated land requirements (demand); (2) available land area (supply) by considering the spatial policies and restrictions; (3) allocation mechanism based on region capabilities.

(1) Demand: the area demanded to cultivate different crops will be determined by the SD sub-model described in the previous section.

(2) Supply: the total available land in each region should be calculated based on spatial policies, restrictions, conversion rules and area of protected lands. For this, we must reduce the total area of restricted lands, protected lands, and non-convertible lands from the total area of agriculturally usable lands.

(3) Allocation mechanism based on capabilities: it is assumed that crop allocation depends on the suitability of each region and current cropping pattern. In fact, crops tend to be cultivated in more suitable areas in long term due to rationality of farmers and ecological enforcements. But these changes could not be happen suddenly because of farmers’ adhesion to current situation and their resistance to change and the need for time to develop infrastructure. Then two factors influence on the future cropping pattern: (1) suitability and (2) change elasticity. In this study, two criteria have been introduced to represent these two factors (see Eq. 1 and Eq. 2) inspired by Balassa index (Balassa, 1965).

RSAi,j determines the relative suitability advantageous of region i for the cultivation of crop j compared to other regions. Similarly, RCE determines the relative change elasticity of region i for the cultivation of crop j compared to other regions. Total relative capability (TRC) of region i for the cultivation of crop j could be calculated by combining these two criteria (see Eq. 3).

Parameter C is the suitability coefficient and denotes relative importance of the suitability factor to the change elasticity factor. The value of C depends on the behavior of people in every society and the time step of simulation. If people’s adherence to the current situation is low and the preferences to move towards optimality are high, then the parameter C will be higher and vice versa. Also if the time step of simulation is longer, then the time needed for change towards suitability is broader, and consequently the parameter C will be greater. This parameter should be calibrated in terms of desirable time step and the conditions of case study before starting the simulation.

In iteration t, the algorithm allocates a percentage of total available land of region i (ali) to crop j cultivation (CLAi,j,t) based on TRC and competitive advantage of crop j (caj,t) (see Eq.4).

Competitive advantage of crop j is “relative advantage of crop j to all other crops”. The caj,t will be determined iteratively. It controls the sum of allocated area to crop j in all regions is equal to the estimated demand from SD model at country level. The value for all crops is 1 in the beginning and iteratively increase for crops that allocated area is smaller than the estimated demand while decrease for crops that allocated area exceeds the demand. The algorithm will be terminated when the total allocated area to each crop is equal to its estimated demand at country level.


Calibration and validation of model

The data of 1996 to 2006 were used to calibrate the SD sub-model in terms of parameters specifying and variable settings. Econometric techniques were used to obtain the response coefficients of SD model (e.g. price response of demand, productivity response to production). In cases where estimation techniques could not be applied (e.g. due to short time series), response coefficients have been embedded from the literature. Initial value of stock variables in SD model will be equal to its value in the first year of simulation.

The main parameter that needs to be calibrated in the crop allocation sub-model is suitability coefficient of each crop. By taking 1998 as the base year, value of parameter C was determined so that the most correlation occurs between actual cropping pattern and simulated cropping pattern in 2006. Results show that the values of C in Iran and for an 8 years’ time step is 0.3, 0.4, 0.5 and 0.3 for wheat, barley, rice and maize, respectively.

There are two ways to validate a SD model: structural validation and behavioral validation (Barlas, 1996). Structural validation is the process of determining whether the theories and assumptions underlying the conceptual model are correct, reasonable, and sufficient. Behavioral validation is the numerical process of determining whether the model behaviors accurately represent the conceptual framework of the model (Qudrat-Ullah, 2012; Mesgari et al., 2017). Since we developed the SD model based on the expert knowledge and the available literature, structural validation of the model has been approved in several focus groups in academic centers related to agricultural and economic studies. Also the model’s behavioral validation is demonstrated by comparing its results to historical trends using data of 2007 to 2014. According to Barlas (1996), a model will be valid if the standard error rate is smaller than 5%. In this research, examined variables include production and cultivated area. Table 2 reports standard errors for the period 2007-2014. Also Fig. 5 shows the dynamic trend of historical values of cultivated area and production variables, compared with their simulated values for wheat and barley between 2007-2014. All the error rates were smaller than 5%, which confirms the validity of the proposed SD model.

Table 2. Average value of the simulation results and historical data of the system dynamics model

Figure 5. Detailed validation results of system dynamics for wheat and barley: a) production, and b) cultivated area.

In addition, the data from2006 to 2014 was used to verify the whole integrated model and to evaluate its ability for projection by comparing its results to historic trends. R-squared values of the regression model between the historical data and simulated data were 0.80, 0.76, 0.61, and 0.76 (p<0.01) for wheat, barley, rice and maize, respectively. These results indicate that the model is reliable and it could be used to project the crop pattern in the future based.

In this study, one probable scenario is defined for the future 25 years in accordance with the combinations of exogenous variables including population growth (medium), economic growth (4.9 %) and precipitation amount (212 mm). This scenario is based on the world population prospects of the United Nations and it almost keeps the present pace with demographic and economic development. Other scenarios could be designed and compared which we leave them to future studies.

Empirical results

Empirical results could be presented in two categories. First, the total cultivated area of each crop that is the output of SD sub model. Second, the spatial distribution of cultivated areas in the geospatial regions. Fig. 6 presents the changes in the area of national cultivated land for wheat, barley, rice, and maize during 2015–2040. It is obvious that although a constant growth in total area of cultivated land is predicted, there is a significant difference in the intensity of variation for different crops. There are almost 25 million hectares of agriculturally usable land in Iran, of which only 16 million hectares are currently under cultivation. Then the predicted increase for cultivated land of different crops will be mainly from two sources: agricultural development in agricultural lands that are not currently under cultivation, and displacement of crop’s farms.

Figure 6. Cultivated area changes at country-level for major crops during the period of 2015-2040

The total area of cultivated land represents a considerable and stable growth for wheat and barley, where the area is predicted to be about 7.7 and 2.2 million hectares in 2040 for wheat and barley, respectively. Most of this increase occurs during the first period of simulation of the model and then, it will continue very slowly. Wheat shows the highest growth rates among the four crops, while the total area of rice and maize will be approximately constant.

Fig. 7 illustrates the national geospatial distribution of cultivated areas for four crops in 2015 and 2040 and in 31 provinces. These maps represent the share of each crop from the total agricultural land area of the province, i.e., the range of 25-50% for wheat means that in the province, 25-50% of the total agricultural land is allocated to wheat. The characteristics of the predicted changes in geospatial distribution of major crops could be interpreted from Fig. 6 and Fig. 7.

Figure 7. Percentage of the agricultural area allocated to main crops in 2015 and simulation results for 2040

Wheat is the most important crop in Iran, traditionally cultivated throughout the country. The growth rate of wheat cultivated lands will be higher in the western mountainous regions and less in the central arid regions. Just like wheat, barley is also a geographically ubiquitous crop; its cultivation area will be the second largest in Iran. Rice is mainly grown in northern Iran, where high temperature and high precipitation are favor to rice cultivation. Of course, some rice can be distinguished in other regions due to historic and the desire of farmers to maintain the previous state. Maize farms have the lowest share of cereal cultivated areas, and the area of its lands will decrease in all provinces. This may be due to the water shortage in Iran in the future.


This paper describes the development of a dynamic model at the provincial and national levels, combined SD sub-model with crop allocation sub-model. This model is able to assess and to analyse the changes of cropping pattern under different scenarios. Validation results demonstrate that the model is reliable for capturing the dynamics of crop demand and cropping pattern. This shows the immense potential of the model to be a decision tool, especially for policy formulations, land use change modeling, and agricultural management. Using the model in Iran’s agricultural system, future demand and spatial distribution of major crops have been calculated up to 2040 with 8th year time step. It can provide theoretical and technical support for agriculture management in Iran.

The suggested model has three major characteristics in comparison with current land use change models:

(1) Using the SD for demand estimation gives a strong to the model in relation to agricultural economics and agricultural system management theories. SD is an excellence way to separate investment and production and demand variables, and to include price effects and trade dynamics.

(2) Embedding SD model in integrated model enables us to predict the possible demand, production, trade, and geographical distribution under different ‘‘what-if’’ scenarios. This advantage allows to rapidly test a wide range of possibilities (Castellazzi et al., 2010). It would be very useful especially when we want to simulate changes of system under the conditions of policy making. Although the base scenario was limited and general without consideration of climate change and detailed government policies, these topics could be considered in the future studies.

(3) Usually in the developing world all the data reported could be fetched at not lower than the district levels and size of those districts also may vary to a greater extent (Priya & Shibasaki, 2001). Since the suggested approach is designed for large provincial scale, it has special advantage for simulation of developing world that the data are always a limitation. However, further researches should be focused on improving the prediction of the model, and the field level interactions within the system.

In summary, this study adds a noticeable contribution to land use change modeling by logically integrating the socioeconomic and spatially explicit data into a SD–CLUE framework. The characteristics of the suggested approach turn it into a useful tool for analysis of agricultural system and policy making. By successful applying the model in Iran, we believe that the model could be used in other countries because Iran is one of the best examples for displaying diversity from one place to another in terms of climate, natural and socio-economic conditions. As a general conclusion, this model not only can be utilized as an effective decision support system (DSS) for the land use planning and management of the agricultural sector, but also can be an appropriate basis to develop rather integrated models for studying agricultural and agronomics systems.


Akıncı H, Özalp AY, Turgut B, 2013. Agricultural land use suitability analysis using GIS and AHP technique. Comput Electron Agric 97: 71-82.

Balassa, B, 1965. Trade liberalisation and "revealed" comparative advantage. The Manchester School 33 (2): 99-123.

Barlas Y, 1996. Formal aspects of model validity and validation in system dynamics. Syst Dynam Rev 12 (3): 183-210.<183::AID-SDR103>3.0.CO;2-4

Britz W, Verburg PH, Leip A, 2011. Modelling of land cover and agricultural change in Europe: Combining the CLUE and CAPRI-Spat approaches. Agric Ecosyst Environ 142 (1): 40-50.

Castellazzi MS, Matthews J, Angevin F, Sausse C, Wood G, Burgess PJ, Brown I, Conrad KF, Perry JN, 2010. Simulation scenarios of spatio-temporal arrangement of crops at the landscape scale. Environ Model Softw 25 (12): 1881-1889.

Doglioni A, Primativo F, Laucelli D, Monno V, Khu ST, Giustolisi O, 2009. An integrated modelling approach for the assessment of land use change effects on wastewater infrastructures. Environ Model Softw 24 (12): 1522-1528.

Forrester JW, 1958. Industrial dynamics: a major breakthrough for decision makers. Harvard Bus Rev 36 (4): 37-66.

Li FJ, Dong SC, Li F, 2012. A system dynamics model for analyzing the eco-agriculture system with policy recommendations. Ecol Model 227: 34-45.

Luo G, Yin C, Chen X, Xu W, Lu L, 2010. Combining system dynamic model and CLUE-S model to improve land use scenario analyses at regional scale: A case study of Sangong watershed in Xinjiang, China. Ecol Complex 7 (2): 198-207.

Mesgari I, Jabalameli MS, Barzinpour F, 2017. System dynamics modeling for national agricultural system with policy recommendations: application to Iran. Pak J Agric Sci 54 (2): 457-466.

Parker P, Letcher R, Jakeman A, Beck MB, Harris G, Argent RM, Hare M, Pahl-wostl C, Voinov A, Janssen M, Sullivan P, 2002. Progress in integrated assessment and modelling. Environ Model Softw 17 (3): 209-217.

Peng M, Chen H, Zhou M, 2014. Modelling and simulating the dynamic environmental factors in post-seismic relief operation. J Simul 8 (2): 164-178.

Pilehforooshha P, Karimi M,Taleai M, 2014. A GIS-based agricultural land-use allocation model coupling increase and decrease in land demand. Agr Syst 130: 116-125.

Priya S, Shibasaki R, 2001. National spatial crop yield simulation using GIS-based crop production model. Ecol Model 136 (2): 113-129.

Qudrat-Ullah H, 2012. On the validation of system dynamics type simulation models. Telecom Sys 51 (2-3): 159-166.

Rosegrant MW, Msangi S, Ringler C, Sulser TB, Zhu T, Cline SA, 2008. International model for policy analysis of agricultural commodities and trade (IMPACT): Model description. International Food Policy Research Institute, Washington DC.

Senge PM, 1994. The fifth discipline: The art & practice of the learning organization. Dell Publishing Group Inc, NY.

Schaldach R, Alcamo J,2006. Coupled simulation of regional land use change and soil carbon sequestration: a case study for the state of Hesse in Germany. Environ Model Softw 21 (10): 1430-1446.

Schreinemachers P, Berger T, 2011. An agent-based simulation model of human–environment interactions in agricultural systems. Environ Model Softw 26 (7): 845-859.

Sharon A, de Weck OL, Dori D, 2011. Project management vs. systems engineering management: A practitioners' view on integrating the project and product domains. Syst Eng 14 (4): 427-440.

Sterman JD, 2000. Business dynamics: systems thinking and modeling for a complex world. No. HD30. 2 S7835. McGraw-Hill Education, Boston.

Verburg, PH, Overmars KP, 2007. Dynamic simulation of land-use change trajectories with the CLUE–s model. In: Modelling land-use change: Progress and Applications, Chapter 18; Koomen E et al. (eds), pp: 321-335. Springer.

Verburg PH, Overmars KP, 2009. Combining top-down and bottom-up dynamics in land use modeling: exploring the future of abandoned farmlands in Europe with the Dyna-CLUE model. Landscape Ecol 24 (9): 1167-1181.

Verburg PH, Soepboer W, Limpiada R, Espaldon MVO, Sharifa M, Veldkamp A, 2002. Land use change modelling at the regional scale: the CLUE-S model. Environ Manage 30 (3): 391-405.

Wang J, Chen J, Ju W, Li M, 2010. IA-SDSS: A GIS-based land use decision support system with consideration of carbon sequestration. Environ Model Softw 25 (4): 539-553.

Zarghami M, Akbariyeh S, 2012. System dynamics modeling for complex urban water systems: Application to the city of Tabriz, Iran. Resour Conserv Recy 60: 99-106.

Zhang Y, Li C, Zhou X, Moore B, 2002. A simulation model linking crop growth and soil biogeochemistry for sustainable agriculture. Ecol Model 151 (1): 75-108.

Zhong TY, Zhang XY, Huang XJ, 2009. Simulation of farmer decision on land use conversions using decision tree method in Jiangsu Province, China. Span J Agric Res 7 (3): 687-698.

Zhou M, Pan Y, Chen Z, Li B, 2016. Enterprise behaviour under Cap-and-Trade conditions: an experimental study with system dynamic models. J Simul 10 (1): 12-23