An Integrated Ecological-Social Simulation Model of Farmer Decisions and Cropping System Performance in the Rolling Pampas (Argentina)

Changes in agricultural systems are a multi-causal process involving climate change, globalization and technological change. These complex interactions regulate the landscape transformation process by imposing landuseandcover change (LUCC)dynamics. Inorder tobetterunderstandand forecast theLUCCprocess we developed a spatially explicit agent-based model in the form of a Cellular Automata: the AgroDEVS model. The model was designed to project viable LUCC dynamics along with their associated economic and environmental changes. AgroDEVS is structured with behavioral rules and functions representing a) crop yields, b) weather conditions, c) economic profits, d) farmer preferences, e) adoption of technology levels and f) natural resource consumption based on embodied energy accounting. Using data from a typical location of the Pampa region (Argentina) for the period 1988-2015, simulation exercises showed that economic goals were achieved, on average, each 6 out of 10 years, but environmental thresholds were only achieved in 1.9 out of 10 years. In a set of 50-years simulations, LUCC patterns converge quickly towards the most profitable crop sequences, with no noticeable trade-o between economic and environmental conditions.


Introduction
. In agroecosystems, land-use and cover change (LUCC) is driven by simultaneous responses to economic opportunities, institutional factors and environmental constraints (Lambin et al. ). Di erent methodological approaches are used for assessing LUCC by capturing the dynamic process influenced by complex interactions between socio-economic drivers and biophysical conditions. For example, the spatially explicit LUCC models assess location suitability for di erent land uses and allocate changes to grid cells based on suitability maps (Verburg et al. a). Another modelling approach has been developed for capturing the behavior of the real actors of land-use change: individuals and/or institutions ("agents") become the objects of analysis, modelling and simulation, paying explicit attention to interactions between these agents of change (Castella & Verburg ). In this modelling approach, farmers' responses to environmental, economic and sociological constraints Region under study and cropping system description . AgroDEVS was applied to simulate LUCC dynamics in Pergamino, a typical location in the Rolling Pampas (Figure ), the most productive subregion of the Pampas where most of the annual cropping is concentrated (Hall et al. ). The Pampa is a fertile plain originally covered by grasslands, which during the s and s was transformed into an agricultural land mosaic by grazing and farming activities (Soriano et al. ). The predominant soils are typical Arguidols and the annual rainfall averaged mm (Moscatelli et al. ). AgroDEVS simulates LUCC using the most frequent crop types in the Pampa region (Manuel-Navarrete et al. ): ( ) the wheat/soybean double-cropping (W/S); ( ) maize cropping (M ), and spring soybean cropping (S). In this paper, the terms "crop"; "crop type" and "land use" (i.e. cropping systems) were used interchangeably as well as "farmer" and "agent". The agronomic decisions (i.e. genotype selection, fertilizer management, pest control, sowing date, and soil type) were representative of the most frequent situation for each of the cropping systems.

Cellular Automata-based modelling and simulation approach .
For the modelling activity we followed a scenario-based analysis, exploring agricultural mosaic dynamics. Then we adopted an agent-based modelling (ABM) approach, by identifying both landscape and agent-specific descriptors as parameters (fixed) or attributes (variable). We defined parameters as any fixed condition for describing the behavior or the condition of a model element. On the other hand, the attributes represent the system changes during the model run period. Lastly, we adopted a formal model-based simulation framework to specify mathematically both the parameters and attributes due to di erent behaviors. For this purpose, we encoded AgroDEVS using the Cell-DEVS formalism. Cell-DEVS is an extension for Cellular Automata of the more generic Discrete EVent System Specification (DEVS) formal modelling and simulation framework. On the one hand, the DEVS formalism permits to express and combine any kind of dynamical system (continuous, discrete event, discrete time) in a mathematical form that is independent of any programming language. On the other hand, Cell-DEVS provides the modeler with a meta-language tailored to facilitate the expressing of systems where the spatial arrangement of "cells", and their behavior, plays a key role. By adopting the DEVS-based approach, AgroDEVS becomes readily linkable with other DEVS models developed by others for di erent domains (e.g. climate, biology, sociology, and economics), potentially using heterogeneous techniques (e.g. di erential equations, equilibrium models, optimization models, stochastic processes). In Figure   .
Behavioral rules are used in the Cell-DEVS language to define the change of attributes that are local to each cell. These variables can express properties for that unit of space (e.g., in the -dimensional case it can be an agricultural plot). When the modeler uses such variables to express also the attributes of "agents" located at cells, then the ABM approach can be merged into the system. In the case of AgroDEVS each agent (namely, a farmer) remains fixed to his or her plot, and rules are used to express changes both for the physical environment and the for farmer, which eventually a ect each other. In this work, we adopted the CD++ simulation so ware toolkit, which is capable of interpreting Cell-DEVS and DEVS models and of simulating them. CD++ implements a generic DEVS "abstract simulator", consisting of a standardized algorithmic recipe that specifies how to simulate any DEVS model independently of any programming language of choice. We believe this generic, reusable and extensible DEVS-based approach provides AgroDEVS with very desirable scalability and sustainability features. The full model represents the collective (i.e., landscape-level) function that emerges from the aggregation of all farmers' outcomes. It also depends on exogenous variables (e.g., climate, crop (output) prices, and production costs) as well as endogenous variables (e.g., the farmer's technological level, the outcomes of neighboring farmers, and each farmer's performance history). The model proved particularly well suited to reproduce empirical situations where (a) there are changes in the relative production/output prices between potential land uses, (b) there is a specific climate regime that impacts on crop yields, or (c) there are varying aspiration levels for the farmers (Gotts & Polhill ). .
In Appendix E we provide more technical details about the model specification approach, its modelling language, rule specification, simulation framework and experimentation environment.

Model verification, validation and sensitivity analysis .
We submitted AgroDEVS to a systematic verification and validation process (Wilensky & Rand ). Firstly, a code walk-through (Stern ) was performed to review the model formulation and to ensure that all design concepts and specifications are correctly reflected in the code. In addition, AgroDEVS was run with very few farmers in the grid ( -) and results were examined closely (e.g., following dynamics of specific farmers, inspecting the outliers). The intrinsic complexities and uncertainties on both the magnitude and the nature of the forces drivingland use change, lead to expanding the scope of the straightforward evaluation between simulated and observed patterns in the model validation phase (Bert et al. ; Le et al. ; Nguyen & de Kok ). Our model development process seeks to validate the ABM through a process of continuous adaptation using feedback from the stakeholders (Ligtenberg et al.
). Specifically, AgroDEVS development entailed a continual discussion process with stakeholders (data not shown) from the study area in order to review both the rules and the assumptions of the model that initially came from the literature review. This approach directly engages stakeholders in model development by embedding it within the social process of policy development (Moss ). This changes the validation problem into an advantage: the agreement of participants or stakeholders may be an indicator of the validity of a simulation model (Troitzsch ). In AgroDEVS, the evaluation of the simulation is guided by the expectations, anticipations, and experience of the community that uses it for practical purposes (Ahrweiler & Gilbert ), and this supports the view that the very meaning of validity is dependent on the purpose of the simulation models under examination (Küppers & Lenhard ). Moreover, it is possible to develop a model that fits the data perfectly with a model structure that does not capture any real processes (Cooley & Solano ). Due to the complex nature of ABMs (e.g., nonlinear responses to parameters) a broad single-parameter SA (Appendix A) was performed (Railsback & Grimm ). This type of analysis assesses the e ect of each parameter over a wide range of values, as opposed to the traditional local analysis, in which each model input is varied by a standard small amount (Ligmann-Zielinska et al. ; Railsback & Grimm ). The goodness-of-fit of each simulated LUCC was calculated using ) Normalized RMSE: v = RM SE/Observed mean %T otal area; ) Probability of a Match: P M = #matches/(#matches + #mismatches); and ) Index of Observed Fit: IOF = (2 × P M ) − 1.
. P M are the chances that a simulation can correctly predict the order of any pair of observations. The number of matches (#matches) is calculated by counting the set of ordered pairs of observations that match the predicted ordered pairs of a simulation (Thorngate & Edmonds ). IOF is an index that derives from the calculations of the PM. Its values range from IOF = +1 (all observations match predictions), through IOF = 0 (half of the observations match), to IOF = −1 (none of the observations match predictions).

Model description
. The model description follows the ODD protocol "ODD" (Overview, Design concepts, and Details) to standardize the published descriptions of individual-based and agent-based models (Grimm et al. ). We describe in the following subsections ) the model's purpose, ) the system variables and ) the process overview to emphasize the main message of the model outputs. The remaining items of the ODD protocol (initialization conditions, submodels description, input data, modelling approach, and design concepts) can be found in the Appendixes B to F.

Model's purpose .
The main modelling purposes (see Edmonds ) of AgroDEVS are ) Prediction, as there is a need to be able to anticipate aspects of the agricultural system that are not currently known; and ) Explanation, as the model tries to identify the causal interactions between the main driving forces of LUCC phenomena (prediction alone does not provide insights on the internal dynamics driving the evolution of system processes).

System variables .
AgroDEVS maps farmers onto a regular grid in order to initialize the simulations. The model consists of two entities: ) the landscape and ) the agents that operate on the landscape. Each entity has its own set of fixed parameters and attributes that evolve throughout simulation cycles. The landscape parameters are a) the number of agents and b) the owner/tenant agent ratio and the landscape attribute is the overall outcome from the integration of all individual agent attributes results within the simulated landscape. The agent parameters are a) the land rental price (RP ) for the tenants, and b) its location on the grid. The attributes for each agent are a) the technological level (T L), b) the crop type allocation (or land use, LU ), c) the economic profit (P ), d) the renewability level (RL) of the embodied energy (emergy) consumption (see Appendix C Renewability level calculations section for an emergy concept explanation), e) the aspiration level (AL), f) the environmental threshold (ET ), and g) the weather growing condition (W GC).

Process overview .
AgroDEVS simulations advance with an annual time step, representing a single cropping cycle (CC) (Figure ).
Crop yields under di erent WGC were previously simulated (Appendix B, Table ) using the DSSAT model (Jones et al. ). We used local management data for defining resource level in each T L (e.g., fertilization regime, genotype) as well as five contrasting historical weather records for defining the di erent W GC levels. attributes and attributes calculated during each CC, respectively. Square boxes denote calculations performed at each stage (for more details of each step see Appendix C). Rectangles represent goal fulfillments, diamonds are conditional rules, and parallelograms are an agent's decisions. Solid and dotted lines denote Control flow and Data flow, respectively. Red borders highlight when the agent performs a neighborhood analysis, and grey background indicates that the calculated value shall be used in the next cropping cycle. .
At the start of the AgroDEVS simulation process (Figure , Stage ), each initial agent's configuration (LU t ) is exposed to a climate-related growing condition level W GC that ranges from very favorable to very unfavorable for high crop yield achievement. Then, based on ) the previously simulated crop yields, ) the crop price, ) the production cost of each T L, ) the land rental price and ) the crop type allocation into the farm household, AgroDEVS calculates the farmer's P t : the profit P for the current time t, or similarly, for the current cropping cycle CC. At this first step, AgroDEVS also calculates RL t (see Appendix C for RL calculation details) as a mea-sure of the renewable energy consumption of the cropping cycle. At the next stage, the initial agent's aspiration level (AL t ) is adjusted by means of the current W GC t level resulting in a new climate-adjusted aspiration level (CAL t ) (Figure , Stage ). Then, AgroDEVS uses P t and RL t values calculated in Stage for assessing the fulfillment of both the environmental (ET ) and economic (AL) goals (Figure , Stage ) of each agent. An agent will keep its crop type allocation once the P t value is greater than or equal to its CAL t at each CC, while nonfulfillment of environmental threshold (ET ) does not alter its crop type allocation decision. AgroDEVS assesses whether the farmer can upgrade or downgrade its T L for the next CC, according to their economic performance in the current CC (Figure , Stage ). For this adjustment, AgroDEVS defines a set of di erent working capital thresholds (W CT ) in order to access di erent T L values. The W CT fulfillment is assessed using P t . If P t is lower than the respective W CT , the farmer must lower its T L (and vice versa). Eventually, agents can remain at the lowest T L, regardless of its P t value (i.e., no agent is forced out of business). During this stage, AgroDEVS also adjusts the farmer's aspiration level (AL) for the next CC. This setting is based on a) the farmer's perception of the W GC level of the next CC, and b) the agent's failure or success at achieving the AL in the previous CC, respectively. Values of ) crop prices, ) adjustment factor of the aspiration level due to W GC, ) adjustment factor of the aspiration level due to T L, and ) working capital threshold for T L are shown in the Appendix B, Tables to , respectively).
The simulated LUCC patterns replicated the overall trend towards soybean-dominated landscapes observed in the region since the mid-s. The ordinal pattern analyses (OPA) showed a similar, relatively high P M for maize, wheat/soybean, and soybean, evidencing the model's capability to predict ordinal (higher/lower) values of crop type cover. However, the accuracy in predicting the magnitude of these changes was lower (Figure ). Based on the RMSE method, the goodness-of-fit of the wheat/soybean simulated cover was the lowest among the three-crop types, resulting in a significant overestimation of this cropping area. The model underestimated simulated soybean cover. In addition, the v values (i.e., the RMSE relative to the observed mean) was remarkably higher for wheat/soybean than in both maize and soybean. Profit, renewability, and technological levels . Profit (P ) and renewability level (RL) values for the Pergamino simulation correlated without any significant trend throughout the study period ( Figure ). Remarkably, in only four years of the simulated period, the simulated landscape exhibited higher RL values than the ecological threshold (ET > 50%RL). Inter-agent variability in the simulated landscape depicted a range of decisions (in terms of crop type allocation) that resulted in a great variability on average P values during the simulation ( Figure   . The interquartile P range ( -% percentile) varied between and US$/ha (Figure , Panel A). However, the RL interquartile range was not as high as the observed in P (Figure , Panel A). The relative homogeneity between each cropping production system, in terms of external inputs, generated a remarkably low variability in RL terms during the whole simulated period. In this case, the interquartile RL range varied from . % to . % with a median value of . %, almost ten percentage points lower than the environmental threshold (ET ) of % (Figure , Panel A). While in terms of RL, the variability between agents was significantly lower than those observed for P , it was possible to detect both maximum and minimum RL values showing the model capability to simulate a wide area of decisions during the simulated period. Minimum average value for an agent outcome in terms of RL was % and the maximum RL average value for an agent was . % (Figure , Panel A), a condition of very close agreement with the environmental threshold (ET ). Agent behavior in AgroDEVS was also assessed by inspecting the intra-agent variability, using both P and RL during the simulation period ( Figure , Panel B). The interannual coe icient of variation of P showed a median value of %, with an interquartile range between . % and . % ( Figure , Panel B). The RL interannual variation was significantly lower than P , exhibiting a median coe icient of variation of . % and an interquartile range between % and % ( Figure , Panel B).

Figure :
Inter-agent and intra-agent variability of the Pergamino -simulation, and the Ecological threshold (ET ) and economic aspiration level (AL) goal agreements. Panel A shows the average profit (P ) and renewability level (RL) of the agents in the Pergamino -simulation. The horizontal solid line shows the overall average of P and RL (n = 625). The extremes of the whiskers represent the % and % quartiles, and the numbers show the minimum and maximum P and RL average values for all agents. Panel B shows the coe icient of variation of Profit (cv P ) and of Renewability Level (cv RL) for all agents. The horizontal solid line shows the overall average of cv P and cv RLV (n = 625), the extremes of the whiskers represent the % and % quartiles, and the numbers denote the minimum and maximum cv P and cv RL values. Panel C shows the average percentage of agreement (% goal agreement) for all agents of the Pergamino simulation. The horizontal solid line shows the overall average of ET and AL goal agreement (n = 625), the extremes of the whiskers represent the % and % quartiles, and the numbers show the minimum and maximum ET and AL goal agreement values. .
An assessment of the simulation results was carried for both the environmental (ET ) and economic (AL) goal agreements exhibited by agents throughout the simulation period ( Figure , Panel C). The goal agreement metric indicates the percentage of years during which an agent fulfills each of the goals (environmental and economic). The agents exhibited an interquartile range of economic (AL) goal agreement between % and % with a median of % and maximum and minimum of . % and . %, respectively (Figure , Panel C). These values are noticeably higher than the RL goal agreement that showed a median of . % with an interquartile range between % and % (Figure , Panel C). The agent's capability for adjusting the AL, based on both the current climate condition (W GC) and the AL fulfillment in the previous agricultural cycle, could explain the better goal agreement when compared against RL, which is a fixed parameter. When inspecting the extreme values for RL goal agreement, it was possible to identify a single agent exhibiting a maximum value of . %. This means that, under the same economic and climate conditions, this agent was able to fulfill nearly out of years the RL threshold (ET ) by means of its crop type allocation decisions.
. Lastly, the model outcome was also assessed by inspecting the final distribution of the three di erent technological levels (T L) across all agents of the Pergamino simulation ( Figure ). The closely related observed and simulated agent distribution underline the model's ability to represent the process of agricultural intensification evidenced by a higher proportion of farmers using a technological level with higher levels of inputs. Both observed and simulated T L distribution patterns exhibit a dominance of the high T L. However, the model retains a greater percentage of agents at the lowest level (T L low) than in the observed data. It is possible that the model structure, which allows agents to remain using the lowest T L despite not reaching the minimum wealth to face those costs, is an explaining factor for the overestimation of Low T L at the end of the simulation. The objective of the long-term simulations was to assess the magnitude of the e ect of the tenure condition (i.e., an owner or tenant-dominated landscape) and climate (i.e., five contrasting climate regimes) on both LUCC (  Profit and Renewability level . The variability between long-term scenarios was lower for both P and RL than for the simulated LUCC ( Figure  ). Regarding the variability induced by the climate regime, the simulated landscapes were stabilized at increasing P -values, as climate scenarios were better (Figure : Panels L, A, and H). This e ect of profit improvement occurred in both tenant-and owner-dominated landscapes. In this latter case, the stabilized P values were due to the di erential income associated with non-payment of land rental. .
Although the simulated landscape configurations were clearly di erent ( Figure ), the RL variability among long-term scenarios under constant climate regimes (when stabilized) showed very small changes (less than %) between the maximum and minimum RL final values (

Discussion
. Determining patterns of land use change is undoubtedly a multifaceted challenge. When relying on simulations, several factors come into play such as climate influence on crop yield, the farmers' decisions, economic prices and costs, and the cognitive description of farmer behaviors (Hare & Deadman ). Thus, the construction of a LUCC simulation model entails necessarily the coupling of social and environmental models (Müller et al. ). In this paper, we developed an ABM expressing several of these characteristics. AgroDEVS simulations were able to reproduce observable LUCC trends of the most representative cropping systems in the region under study. An ABM validation process has the peculiarity of being subjected to conflicts between achieving accuracy in matching the outcome of a simulation or in the processes simulated (Brown et al. ). This trade-o is usually solved based on the research goals. In the case of LUCC simulation models, both aspects are important. Predicting the LUCC trends is extremely important in decision making by policy makers (Verburg et al. b). However, this information should be supplemented through an understanding of the underlying LUCC processes involved. This is required to identify potentially unsustainable land use regimes and correct them. AgroDEVS' results show that its structure is able to detect the overall trend on land use changes. This is done through clear and explicit modelling that reflects key process dynamics such as the climate influence on crop yields, farmer decisions and the landscape emergent properties due to farmer interactions at smaller scales. .
The contribution of the model developed in this work can be assessed by analyzing carefully the LUCC simulation results for Pergamino -, while considering the trade-o s between output accuracy and processes understanding. Although the model was able to simulate the land use change dynamics of the three crops analyzed (i.e., the ordinal fit is always higher than . ), the adjustment based on the distance between observed and predicted (i.e., the v value) could be improved by the inclusion of other variables or exploratory processes (e.g., agricultural policy decisions not considered; the dynamics of prices and costs, etc.). However, the current model structure maintains the relative profit between activities which is highly sensitive to environmental conditions (W GC). Thus, the distance between the observed and predicted LUCC can be used as a predictor of the di erence between what could have been a LUCC trend (based solely on the response to the environment in order to maximize profit, i.e., LUCC simulated) and another path that did not strictly follow the parameters of maximizing profit (i.e., LUCC observed). Notably, the region studied has been frequently subjected to decisions in agricultural policy (e.g., di icult marketing of some crops, imposition of export taxes) that strongly influence the changes in land use through mechanisms not directly related to profit activity (Porto & Lodola ).
. From the formal modelling and simulation point of view, the DEVS approach (and its related Cell-DEVS spatially specific flavor) o ered several salient features, already discussed in the model description section. Notably, in the context of socio-environmental systems, there is a feature that stands out. Namely, the input-output portbased hierarchical modularity permits the design of interdisciplinary models by composing complex systems through the interconnections of simpler ones. For instance, in AgroDEVS, the DEVS Atomic Model representing climate dynamics can be replaced by another, more accurate or sophisticated DEVS model developed by experts in the climate domain, without requiring to alter the landscape (Cell-DEVS) portion of the system. This approach is consistent with a current trend towards Systems of Systems-oriented modelling and simulation (Zeigler et al. ). This is especially relevant for socio-natural-economic systems, whose domain-specific submodels are constantly subject to revisions, improvements or replacements. Regarding the use of Cell-DEVS to spatially represent farms, abstract cellular landscapes have already been used in other ABM-LUCC models (e.g. FEARLUS, AgroPoliS, Pampas Model; Polhill et al. ; Happe et al. ; Bert et al. ). Results from AgroDEVS indicate that this approach has proven to be an e ective and simple way of spatially modeling agricultural systems. .
Incorporating environmental assessment in LUCC simulation models is a very desirable feature that is beginning to be explored (Veldkamp & Verburg ). The analysis of environmental impacts on managed ecosystems has o en been applied based reductionist approaches, identifying changes very accurately, but reducing the relevance of the results by not addressing an integrated or holistic approach for ecosystems modelling (Shanmuganathan et al.
). The AgroDEVS structure acknowledges the need for a systemic modelling by including an emergy (embodied energy) renewability level as a proxy for assessing ecosystem functioning of the cropping systems studied, an approach that has been tested before with agricultural systems both in the studied region and in other agricultural ecosystems (Agostinho et al. ; Dewulf et al. ; Ferraro & Benzi ). Concerning the use of environmental work, the simulated ecosystems did not show a clear dynamic of increase or decrease in its reliance on energy from the economic system. Moreover, the results from the simulations, both in the Pergamino -and the long-term runs, also failed to show a clear trade-o between the environmental performance (assessed through renewability level) and economic performance (assessed through economic profit) in the studied systems. The simulated dynamics showed that the systems reduced their dependency on exogenous emergy (i.e., fossil energy) to the extent that environmental conditions improved, and thus natural resources become more important. Clearly, this happens because at increasing levels of technology, increased use of external emergy is roughly proportional to the increase in capturing endogenous emergy associated with growth conditions (W GC) that improve crop yields. The high emergy return seems a characteristic of the farming systems of the Pampas region at the field scale (Ferraro & Benzi , ); but the AgroDEVS simulations showed that this pattern is maintained at the landscape scale due to the emergent properties that arise from the integration of individual farmers behavior. .
In the context of changes in agricultural policy in the area under study, the results of AgroDEVS seem to indicate the need for policy options that alter the relative prices of crops, to avoid the predominance of monocultures or systems highly dependent on the soybean crop. In this context, the AgroDEVS capability for identifying potential trade-o s between di erent agroecosystem domains (i.e., economic, environmental) is extremely important in the diagnosis of agricultural sustainability (Tittonell ). Policy options, which have direct e ects on model variables such as prices, can be readily tested into AgroDEVS for their repercussions on farm incomes and LUCC emergent patterns. Moreover, there is always the chance to include other variables in the modelling exercise. In the interest of better representing the heterogeneity amongst di erent farmer's decision logics in AgroDEVS, future e orts could be aimed at exploring the farmers' decision-making process. This exercise could expose new relevant variables that improve the representation of agents' behavior. However, this inclusion would require a new setting and a more complex numerical validation. The cost-benefit balance of these additions should be analyzed carefully to avoid incurring in an overfitting, conspiring against the understanding of the true underlying phenomena under study (Brown et al. ).

Conclusion
. The model presented in this work represents an e ort to integrate, within a complex system simulation, the e ects of weather on crops, the farmer decision logic, and the cropping system profit as main LUCC drivers. AgroDEVS simulations based on real landscape data showed that LUCC direction was better represented than its magnitude, in terms of the land cover dynamics. The long-term simulations showed a dominance of cropping systems that included soybean crop, and this dominance was stronger for monospecific soybean crop in scenarios under constant climate. The double cropping W/S dominated mainly in scenarios under a variable climate. When assessing the farmer condition e ect (i.e., tenure and climate) on LUCC, tenure resulted in much less e ect on LUCC than the weather conditions (W GC). Finally, simulations showed no trade-o s between environmental and economic outcome both in simulation used to validate the model and long-term scenarios. The results suggest that LUCC modelling and its environmental and economic consequences is feasible and useful using an ABM approach. The AgroDEVS simulations allow not only to speculate on the LUCC dynamics, but to gain a greater understanding of the underlying processes involved. Future research should focus on improving the model structure to include di erent agent behaviors (e.g., multiple agent's profiles) as well as social and political factors both for predicting the LUCC direction and to assess their relative magnitudes more accurately. Moreover, incorporating di erent productive regions of Argentina into AgroDEVS could expand the conclusions achieved in this paper, revealing di erences in LUCC dynamics at an ecorregional level.
ulator is CD++, a general purpose discrete-event simulator written in C++ by several generations of developers. The version recommended for this work is: https://github.com/SimulationEverywhere/CDPP_ExtendedState s-codename-Santi.To access the model please contact the authors.

Appendix A. Sensitivity Analysis
In a long-term simulation ( years), we tested the e ect of ) crop prices, ) W GC, ) the Owner/Tenant agent ratio and ) rental price. We run AgroDEVS varying one parameter at a time over a range of values (Table ), while keeping the other parameters at their reference value (i.e. the initialization conditions of the Pergamino simulation; for more details see Appendix B). The range of values for the crop prices was extracted from historical data and are the lowest and highest price of each crop for the years -.

Sensitivity analysis
The WGC scenarios were created by alternating % of the campaigns with historical values (i.e. actual weather for that year) and the other % with the WGC level desired. Thus, five scenarios were created comprising % of historical values and % W GC (V U , U , R, F and V F ). The owner/tenant agent relation ranged from % of the agents being owners and % tenants to % owners and % tenants in increments of % (i.e. / , / , / , / , / ). The rental price ranged from US$/ha , to US$/ha , (Table ).
The sensitivity analysis showed that the model is sensitive to wheat and soybean price changes in terms of profit. When the price of one of these crops rises, the land use cover that dominates the grid is W/S. This crop sequence has a high profit potential, and thus explains the change in profit observed. However, if the soybean price drops below the reference value the farmers cannot compensate with another LU and the profit drops. In terms of renewability, the model was not sensitive to changes in crop prices and all crops behaved similarly. Only a large increase in wheat price caused a % reduction in the RL (Figure ).  The model's response to di erent weather patterns was clear in terms of change in profit (Figure , le ). Under favorable W GC, the profit was higher than the reference value. In the best-case scenario ( % of the years under V F weather), the profit was increased by % of the reference value. Conversely, under unfavorable W GC the profit tended to drop. In the worst-case scenario, ( % of the years under V U weather) the profit decreased %. This behavior is explained by the quick dominance of a more profitable LU (W/S) and management level (High) under favorable weather conditions, and the dominance of a defensive LU such as soybean under unfavorable weather conditions. The model showed little response in terms of RL under di erent W GC (Figure ,  right). The results indicate that under favorable weather conditions the RL could drop slightly due to the dominance of W/S. On the other hand, unfavorable weather conditions would lead to an increase in S, provoking an increase in RL. The model was slightly sensitive to changes in rental cost, leading to profit variation (Figure , le ). As expected, an increase in rental cost diminished the profit in tenants resulting in a decrease in the average landscape profit (no dominance of LU was observed here). On the other hand, the RL was not sensitive to changes in rental cost (Figure , right). This could be explained by the relatively low fraction of the profit that the rental cost represents.

Appendix B. AgroDEVS Initialization Conditions
Crop yields are based on simulations using crop models in the Decision Support System for Agrotechnology Transfer (DSSAT) package (Jones et al. ) that have been calibrated for the studied location (Mercau et al. ). Crop management of each T L (e.g. genotype, sowing date, fertilizer rate, pesticide use, tillage operations, and production costs) was used for running the DDSAT crop yield simulations, and was based on local management data. Crop yield variability during the DSSAT simulated period ( -) was used to obtain the crop yield under contrasting W GC levels. The five W GC levels (very unfavorable, unfavorable, average, favorable, and very favorable) correspond to di erent percentiles (i.e. , , , , and , respectively) of yields simulated using historical weather records. DSSAT simulations do not account for important factors such JASSS, ( ) , http://jasss.soc.surrey.ac.uk/ / / .html Doi: . /jasss. as weeds, diseases, and pests. Thus, we empirically adjusted the attainable crop yield (van Ittersum & Rabbinge ) resulting from DSSAT simulations in order to model actual crop yield. Local data of simulated versus observed crop yield were used for obtaining the adjusting coe icients (attainable to actual yield) at each T L for each crop species (Mercau et al. , ; Satorre et al. ). AgroDEVS reflects three incremental T L: low (L), average (A) and high (H) technological levels; and three crop types (maize, soybean and wheat/soybean double-cropping, see region under study and cropping system description section for production system description). Using DSSAT and historical records of production costs, two of the look-up tables are built: ) a crop yield matrix (Table ) and ) a cost matrix for representing the full crop type ( levels), T L ( levels) and W GC ( levels) combination (Table ). As for crop yields and production costs, the emergy accounting method is used for building a third look-up table for representing the full combination of crop types ( levels), T L ( levels) and W GC ( levels) ( Table ). Historical prices of maize, soybeans, and wheat, as well as production costs (i.e. fertilizers, seeds, pesticides, and harvest and sale costs), were extracted from the Argentine trade magazine "Márgenes Agropecuarios" (http://www.margenes.com). In all scenarios (i.e. the Pergamino simulation and the long-term simulations) we assumed constant output prices equal to median for -.

Maize Soybean
Wheat/ Soybean Table : Adjustment factor αAL(W GC) of the aspiration level (AL) due to weather growing condition (W GC) level. W GC classes: Very Unfavorable; Unfavorable; Average; Favorable; Very Favorable.

Very Unfavorable
Unfavorable Average Favorable Very Favorable -.
-. . . . Table : Adjustment factor αAL(BN ) of the aspiration level (AL) due to agent technological level (T L) and the best neighbor technological level T L t (BN). A positive (or negative) αAL(BN ) value indicates an increase (or decrease) in AL since the best neighbor exhibits a higher (or lower) T L than that evaluated by the agent. -.
-. . Table : Working capital threshold (W CT ) for each technological level (T L) class (expressed in US$/ha.). W CT values equal to % of the highest production cost for each T L, considering and average indebtedness rate of % (AACREA ). T L classes: H (High); (A) Average; L (Low).

Appendix C. Submodels Description
Profit calculations P is calculated as the gross income (yield times product price) minus direct costs. Direct production costs include fixed and variable components. Fixed direct costs do not depend on crop yield (e.g., seed and agrochemicals). Oppositely, variable direct costs are yield-dependent (e.g., harvest, marketing fees and grain transportation). Finally, each farmer P is calculated using the individual crop profit a ected by each crop type allocation into the agent farm. Figure is an excerpt of the rule that calculates P : an agent has a defined LU and T L and is capable of sensing the WGC, and therefore it can calculate its P by multiplying its crop type allocation, yields and prices, and subtracting costs. Figure : Excerpt of the rule that performs agents' profit calculations. In the cases of yield, price and cost, #macro involves lookup tables.

Renewability level calculations
AgroDEVS calculates Renewability Level (RL) values by using the emergy synthesis procedure (Odum ). Briefly, the emergy accounting methodology tries to account for both the natural and human-made capital storages. The emergy accounting method values these storages using a common unit of reference, the solar equivalent joule (seJ). The method accounts for the environmental support provided directly and indirectly by nature to resource generation and processing; it focuses on valuation of the intrinsic properties of ecosystems (Mellino et al. ). For further details on emergy synthesis methodology see Brown et al. ( ) and Ferraro & Benzi ( ). AgroDEVS uses the renewability level as a sustainability metric (Giannetti et al. ). RL is then calculated as the ratio of renewable emergy to total emergy use, as it follows: where, RL(%): renewability level value, R(seJ/ha −1 year −1 ): renewable flows from nature, N (seJ/ha −1 year −1 ): nonrenewable flows from nature, F (seJ/ha −1 year −1 ): imported economic flows, S(seJ/ha −1 year −1 ): services.
Therefore, at each model time step (i.e. a CC) each individual farmer renewability level (RL) is used as a measure of environmental impact. In the figure below ( Figure ), an agent uses its LU and T L along with the sensed WGC to obtain the final farmer's RL value (by a ecting the individual crop renewability level with the crop type allocation into the agent's farm).

Figure :
Excerpt of the rule that performs agents' renewability calculations. In this case, the #macro(emergy) involves emergy lookup tables.

Agent decisions
Both P and RL values for each agent are related to the AL and ET thresholds, in order to trigger the agent's decisions. In AgroDEVS, the economic goal (AL) is dynamic and it is based on the aspiration level adjustment. It represents the currently dominating economic paradigm, where the receiver (the agent) is the market actor who decides the system outcome value (Grönlund et al. ). Oppositely, the environmental threshold (ET ) is fixed and it represents a strong sustainability view, where the value approach is grounded in systems science rather than economic science, where a value focused on the system level is accepted (Grönlund et al. ). During the model simulation process, the fulfillment of the economic goal drives the crop type allocation. Thus, the LUCC process is triggered when the economic threshold (i.e. aspiration level) is not accomplished.

Aspiration level adjustment
A first AL adjustment is based on the W GC, as it follows: where CAL t = climate-adjusted aspiration level (US$/ha), AL t = aspiration level (US$/ha) calculated at the end of the previous CC, αAL(W GC) = adjustment factor of AL due to W GC level for the current CC (see Appendix B, Table for αAL(W GC) values used in the simulations).
A second AL adjustment defines the aspiration level for the next CC and it is based on the learning and adaptation model (Bert et al. ). If P t > CAL t then the next aspiration level is subjected to an incremental adjustment using a weighted average, as it follows: where AL t+1 = aspiration level (US$/ha) for the next CC, CAL t = aspiration level (US$/ha) a er current CC adjustment, P t = Profit (US$/ha) calculated in profit calculations section. Figure : Excerpt of an agent's aspiration level adjustment. In this example, the agents' profit was larger than its AL, so it proceeds to perform an incremental adjustment.
In the figure above we can see an excerpt of the rule for adjusting incrementally the aspiration level using a weighted average when an agent's P t > CAL t . When both the agent and its neighbors finish with the calculation of their P , RL and the fulfillment of both economic and environmental goals, the agent proceeds to adjust its AL for the next cropping cycle. When the economic outcome is lower than the aspiration level (i.e. P t < CAL t ), the agent's perception is extended to include the influence of the physical (Moore) neighbors, and the AL t+1 is adopted by inspecting the neighbors' P outcomes. Thus, the farmer adopts AL t+1 using the neighbors' P data. If there is at least one neighbor under the condition of P t > CAL t then the farmer selects its best neighbor (BN ), in profit terms, and the AL t+1 is calculated using both the AL for the best neighbor, CAL t (BN ), as well as an adjustment factor due to di erences in T L between farmers: where: AL t+1 = aspiration level (US$/ha) for the next CC, BN = the agent with the highest P value in the Moore neighborhood (n = 8), CAL t (BN ) = aspiration level (US$/ha) for the BN in the current CC, αAL(BN ) = adjustment factor of AL due to BN technological level (see Appendix B, Table for αAL(W GC) values used in the simulations).
In the case that P BN ≥ CAL t (i.e. no neighbor meets the economic threshold) then the next aspiration level is subjected to a detrimental adjustment using a weighted average, as it follows: where: AL t+1 = aspiration level (US$/ha) for the next CC, CAL t = aspiration level (US$/ha) a er current CC adjustment, P t = Profit (US$/ha) calculated in Profit calculations.
Di erent values for incremental and detrimental adjustment ( . and . respectively) were applied in order to simulate the farmer willing to tolerate higher payo s more rapidly than lower ones, thus showing greater resistance to downward changes (Gilboa & Schmeidler ).

Technological level adjustment
Farmers are also able to upgrade or downgrade their technological level (T L) a er inspecting their own P outcome in the CC. The rules for adjusting the T L are based on the working capital threshold (W CT ) which represents the highest production cost within each T L ( Table ??). Rules for T L adjustment are: where, T L t+1 = Technological level for the next CC (H: high; A: average and L: low), W CT (T L i ): working capital threshold (average or high) for T L i (low, average and high).
In Figure we can see an excerpt of the rule for adjusting an agent's T L. When both the agent and its set of neighbors are done with calculating their P , RL and evaluating the fulfillment of both economic and environmental goals, the agent proceeds to adjust its T L for the next cropping cycle by comparing its P with the W CT . Figure : Excerpt of an agent's technological level adjustment. In this case, #macro (wc max ) involves WCT lookup table (Appendix B, Table ).

Land use and cover change (LUCC)
Final farmer decision defines the land use configuration. The LUCC process is triggered when the economic threshold is below the economic outcome and also there is at least one neighbor in the Moore neighborhood that meets the economic threshold as follows: where, P t = Profit (US$/ha) calculated as in profit calculations section, P t (BN ): Best Neighbor Profit (US$/ha), CAL t = Aspiration level (US$/ha) a er current CC adjustment.
A er inspecting this condition, the agent selects the crop type allocation of the best neighbor (BN ) as follows:  In the figure above there is an excerpt of the rule that copies the best neighbor's land use allocation. If both an agent and its neighbors have already calculated their P and RL then the model will proceed to select the JASSS, ( ) , http://jasss.soc.surrey.ac.uk/ / / .html Doi: . /jasss. neighbor with the highest P and copy its LU allocation, AL and T L to be used later in the simulation and, if necessary, assigned to the agent.

Pergamino -simulation
Input data from -were used for both the model initialization and to test the model outcomes against the actual land use changes in the studied area. Table presents a summary of these initial conditions.
The observed LUCC dynamics were obtained from agricultural surveys (SAGPyA ), and both the initial ( ) and final ( ) T L frequency distribution among farmers were obtained from the national agricultural census (BOLCER ; INDEC ). In order to represent the main tenure regimes in the studied area (owners and tenants), AgroDEVS was initialized using the Owner/Tenant relation (INDEC ). The land rental price was set to historical values of . t of soybean per hectare (Margenes Agropecuarios ). Crop yield and renewability level values under di erent weather conditions (W GC) and technological levels (T L) were calculated as it is explained in Appendix C profit calculations section. In the initial landscape configuration, each farmer was assigned an initial working capital (W C) according to his initial T L. Initial AL level were fixed for each farmer in order to account for the % of the direct costs (i.e. 0.6W CT ) of the T L adopted in the first CC of the simulation (AACREA ). Oppositely, the ET value remains constant during the simulations and is fixed in RL = 50% (Table ). The accuracy of the simulations was assessed using both ) the squared distances between a simulation's outputs and a set of observations (RMSE) and ) the ordinal pattern analyses (OPA) (Thorngate & Edmonds ). OPA indicates the topological fit between observed and simulated outputs but does not consider their closeness.

Long-term scenarios
AgroDEVS was also run over a -year period, under contrasting scenarios due both ) five W GC regimes (constant unfavorable, constant average, constant favorable, a see-saw pattern of very unfavorable-average-very favorable, and a random regime), and ) two tenure regimes based on the landscape pattern of Owner/Tenant agent relationship ( / and / ). As tenants are more focused on short-term income and are less likely to invest in longer-term management strategies than owners (Soule et al. ), we used AgroDEVS for testing the hypothesis that tenant farmers are less likely than owner-operators to adopt crop allocation decisions that lead to sustainable LUCC trajectories. All crop types and T L were set to equal distribution in the landscape (i.e. % of the total area for each crop and each T L) at the initialization. However, the internal assignment of each crop type for each agent was set randomly. Model simulations were inspected in terms of the dynamics of ) crop type allocation of total area, ) profit and ) renewability level.

Appendix E. Modelling Approach
In Figures and we provide an illustrative, non-exhaustive sample view of the Cell-DEVS specification used for AgroDEVS. Most of the information declaring the structure, components and behavior is found in a text file with .ma extension. Figure : Sample excerpt of an AgroDEVS.ma definition file (Le ) Schematic of the AgroDEVS system. The "climate" DEVS model is connected with the "landscape" Cell-DEVS model (Right).
In Figure we can see the main statements defining structure, components and interconnections through input-output ports. In Figure we provide an excerpt of the agrodevs.ma file highlighting relevant lines. Note that three dots ". . . " in the figures denote lines omitted for the sake of brevity. In Figure  we can see main statements defining behavior for the landscape cellular automaton. Within each cell, a list of rules is evaluated sequentially in a top-down fashion, using a Value, Delay, Condition structure: the first rule that evaluates its Condition to true, will apply the Value to its attributes a er a Delay amount of time. Once a rule is applied, the simulator recommences the cycle evaluating the rules from the top. This mechanism is applied asynchronously, simultaneously and in parallel to all cells in the model. The global timing for the whole cellular space is an emergent property driven by the local timing applied by each independent cell. The pairs (X, Y ) denote the positions for neighbors relative to each currently evaluated cell denoted with (0, 0).

Simulation framework and so ware experimentation environment
As discussed previously in the modelling approach section, AgroDEVS uses a DEVS-based formal approach. One salient feature of the DEVS formalism is the strict and clear separation between simulation algorithms and model specification. As a specialization of DEVS for Cellular Automata, the Cell-DEVS formalism inherits this separation. Di erent Cell-DEVS-capable simulators should be able to simulate a given Cell-DEVS model. In turn, di erent interactive user interfaces can be used to help with the design and maintenance of model specifications and to interact with the simulation exercises. Figure   The CD++ generic simulation engine for DEVS and Cell-DEVS models resides in a Simulation Server. It can retrieve and simulate Cell-DEVS formal specifications of models that are stored in a Database Server. Such models can be written directly by specialists in the Cell-DEVS formalism or by higher-level experts in the agricultural domain using a friendly web-based interface. Once models are defined, they can be reused to launch as many simulation experiments as required, by invoking the simulator also from a web-based experiment interface. Figure describes sequences of steps for two typical use cases: Model Specification and Simulation, denoting the components involved in each of them.

Objectives
Agents pursue to achieve a P level in each CC above their AL while satisfying (or not) a fixed RL. If P is below agents' AL, they will be unsatisfied and will seek a di erent crop type allocation in the agents' neighborhood. If an agent's capital availability drops below the W CT for each T L, they must adopt a lower-cost T L, but agents never quit farming. Both landscape and individual RL emerge from the crop allocation decided for each agent, but it is not used as a farmer s goal. Rather, this is an emergent property due to local rules.

Prediction
Agents predict the future consequences of their decisions (i.e., they build expectations) about P based on past outcomes and current weather information (W GC), following the theory of adaptive expectations (Shell & Stiglitz ). Agents revise their AL in each CC using adaptive rules based on a) the dissimilarity -in the previous CC-between the predicted and the observed P , and b) the seasonal climate forecasts (W GC) for the current CC (see sub-model section for algorithms adjusting AL).

Sensing
Agents have information about their capital availability and consider this variable in their decisions about the potential adoption of higher (or lower) T L. In addition, they are aware of the P achieved by their eight neighbors (Moore neighborhood) during each CC. Finally, farmers are informed about the expected status of external contextual factors (i.e. seasonal forecast) for the current CC.

Agent-Agent interaction
Agent interaction is based on the imitation of the crop type allocation of the best neighboring agent (i.e., the highest P achieved in the neighborhood) in the case of dissatisfaction with its own P achieved (as compared to its AL). Otherwise, the agent repeats its current crop type allocation in the next CC.

Agent-environment interaction
Crop yield simulations reflect the interaction between farm decisions and both climate and soil conditions. These models capture the T L e ect (through increasing fertilization and pesticide use) as well as the W GC of each CC (mainly global radiation and precipitation) and land quality factors (which are constant during the simulation). As crop simulation models represent biotic-unconstrained yield potential, final yields are adjusted in each combination of crop type-T L using potential weeds and pest infestation losses (see simulations input section for an example of crop yield data). Agent decisions have no e ect on resource dynamics, through crop type allocation, and it is only on the current W GC that agents receive feedback about changes in resource conditions.

Stochasticity
Stochasticity is used to assign the T L and the crop type allocation to each agent only at the initialization procedure.

Collectives
Agents do not form or belong to aggregations that a ect or are a ected by other agents in the model version described here. Further model refinements should include collectives.