An Agent-Based Crop Model Framework for Heterogeneous Soils

Climate change and the efficient use of freshwater for irrigation pose a challenge for sustainable agriculture. Traditionally, the prediction of agricultural production is carried out through crop-growth models and historical records of the climatic variables. However, one of the main flaws of these models is that they do not consider the variability of the soil throughout the cultivation area. In addition, with the availability of new information sources (i.e., aerial or satellite images) and low-cost meteorological stations, it is convenient that the models incorporate prediction capabilities to enhance the representation of production scenarios. In this work, an agent-based model (ABM) that considers the soil heterogeneity and water exchanges is proposed. Soil heterogeneity is associated to the combination of individual behaviours of uniform portions of land (agents), while water fluxes are related to the topography. Each agent is characterized by an individual dynamic model, which describes the local crop growth. Moreover, this model considers positive and negative effects of water level, i.e., drought and waterlogging, on the biomass production. The development of the global ABM is oriented to the future use of control strategies and optimal irrigation policies. The model is built bottom-up starting with the definition of agents, and the Python environment Mesa is chosen for the implementation. The validation is carried out using three topographic scenarios in Colombia. Results of potential production cases are discussed, and some practical recommendations on the implementation are presented.


Introduction
Increasing food production involves efficient use of resources including arable land [1]. However, not all the land available for agricultural activities is found in locations with smooth topography and freshwater availability for irrigation. This poses two challenges: the use of the land and the rational use of water. Land irrigation is a major concern, as it represents the largest consumption of freshwater all over the world [2]. Several attempts to improve the irrigation schemes have been proposed over the years, but a systemic understanding of water effects on crops is necessary and involves the use of software tools to achieve this purpose. As it takes a full season to assess an irrigation policy, research has been essentially carried out through numerical simulation. On the other hand, the land available for agriculture (e.g., wasteland), lacks useful information to forecast a food production scenario. In developing countries there are no detailed soil studies available. Then, a first approach is to assume the homogeneity of land based on a limited sample of the soil and to use traditional crop-growth models. Therefore, the crop-growth models assume homogeneous blocks of soil, which makes sense for small regions but is unrealistic at wide regional scale [3][4][5]. Improving these models is not straightforward due to the diverse methodologies and assumptions used to build them. This difficulty was evidenced in the first work of the Agricultural Model Intercomparison and Improvement Project (AgMIP), a global initiative aimed at comparing 27 main models and looking for possible multi-model ensembles [6]. The results reported by Martre et al. [7] indicate that heterogeneity of land at regional scale is a current issue in agriculture.
Some crop simulation models consider soil heterogeneity only with respect to depth, by assuming that the dynamics of the movement of water in the soil only depends on the thickness of the different layers (e.g., AquaCrop [8], STICS [9] and APSIM [10]). However, soil variability across the surface has not been addressed by these models. To face these variations, this work proposes an agent-based model (ABM), that considers a mechanistic approach to crop growth and soil water mass balance. The model proposes a reduced set of equations for an agent to be scalable, and includes uncertainty by adding Gaussian noise in the states and the measurements. The model is fed with climatic variables, irrigation inputs and data from measurements of state variables. With these elements, this model starts from the global behaviour of the crop and allows interpreting the behaviour of smaller portions of land and its response to environmental variables and the amount of water in soil, therefore to propose or to assess irrigation policies.
The application of ABM is not new in agriculture. For instance, an agent structure has been proposed to represent the management decisions of farmers [11,12]. In relation to soil description, [13] presents an ABM framework to model the movement of water through a soil column. To the best of our knowledge, however, no integrated model has addressed surface land variability yet. Instead, large scale models are run in independent simulations and results are collected by a decision support software, e.g., DSSAT in [14]. The topographic surface land variability comes with additional water dynamics that are underestimated by the aforementioned models. To address this issue, a coupling between hydrology and crop growth models was reviewed by Siad et al. [15]. The main conclusions of their work are that coupling models are still at an early stage because data is scarce, and the assessment outputs are reported as regional average.
The main contribution of this work is to propose an agent-based model to represent the complex dynamical interactions between irrigation policy, heterogeneous soil, and crop growth. Each agent represents a portion of land and is associated with a dynamic model of crop growth. The model considers the negative effect of excess or lack of water through an original stress index. Another agent describes the irrigation system, and various irrigation strategies can be assessed in simulation taking constraints on water supply into account. The advantage of the ABM formulation is that it allows flexibility in the definition of the land topography and the consideration of uncertain parameters. In view of the limited (ground or aerial) instrumentation, the control structure involves an estimation scheme based on extended Kalman filtering. The irrigation strategy is based on the optimization of a performance index related to biomass production under irrigation constraints. The software application is illustrated with three topographic scenarios in Colombia. Indeed, 35% of the Colombian territory is available for agriculture but only 7% is currently cultivated. More than 50% of the cultivable area is located on mountainsides, in valleys and on uneven grounds (non-smooth topography and waterlogging in 42%, and 20% of cropping lands, respectively). Moreover, there are sustainable freshwater sources, but they need to be exploited rationally.
This paper is organized as follows. The next section introduces the agent-based model structure, including the estimation scheme and optimal irrigation policy, as well as the implementation in the Python environment Mesa. Three land topologies and four datasets corresponding to different rain patterns are introduced in Section 3, which are used to illustrate the predictive capability of the model and the effect of irrigation on soil stress and biomass production. Conclusions and perspectives are drawn in Section 4.

Model Structure and Implementation
Agent-based modelling (ABM) is a methodology used to simulate interactions between autonomous individuals [16]. An agent is a self-contained, self-directed, modular, and uniquely identifiable individual. The modularity requirement implies that an agent has a boundary. Moreover, every agent is an entity located in a specific environment and capable of autonomous actions, in order to meet some objectives. The complete definition of an agent involves a set of states that varies over time, a set of non-uniform attributes, asynchronous interactions, and uncertainty related to parameters or states [17]. Furthermore, an agent has at least two levels of interaction: an upper layer to interact with other agents, and a lower layer to solve internal processes. In the first layer, the goal is to interact with other agents, and the second deals with the physical world.
The modelling process can be top-down or bottom-up. The bottom-up strategy is selected in this study because one of the modelling goals is to understand how the individual components affect the all-in-all system behaviour. The emergent system behaviour is the result of agent interactions [18]. A general structure (shown in Figure 1) of an agent-based model has typically three elements: 1. A set of agents, their attributes, and behaviours (lower layer). 2. A set of agent relationships and methods of interaction (i.e., interaction topology or upper layer). 3. The agent environment (i.e., spatial topology). For the particular case of this application, two types of agents are considered: (i) portions of homogeneous soil with a surface of regular shape, and (ii) irrigation equipment. The soil agents are associated with a dynamic model representing local crop growth as a function of water availability. The decision or upper level of all agents ensures the respect of the conditions for the exchange of water between neighbouring agents. Crop-soil agents require irrigation at regular time intervals and the irrigation agent answers this demand under the constraints of an irrigation quota.
The global model is intended to predict the evolution of the water content in soil using local information and assess the biomass produced by applying different irrigation schemes. The model structure is shown in Figure 2.
The soil-crop agents are described by a set of state variables x n (k), which are influenced by the climatic variables w(k) and by the irrigation inputs u n (k). Part of the state variables can be measured by sensors with output z n (k) (ground sensors, aerial or satellite measurements) under the influence of measurement errors (or noise) v n (k). The irrigation agent acts as a controller to determine the water delivered to the crop-soil agents. The control action is determined based on information that it receives at regular time intervals k regarding the amount of water available for irrigation of the entire field r(k). As instrumentation is limited, an estimator is required which provides state estimatesx n (k). - Block diagram of the model with n = 1, ..., N crop-soil agents. All variables are described in Table A1 of Appendix A. Each of the model blocks are described in the following sections and the compendium of the variables is presented in Table A1. Subsequently, the Python programming environment Mesa is described, as well as the methodology to define agents in their spatial and temporal environments.

Crop-Soil Agents
The lower level behaviour of every crop-soil agent is related to its internal mechanistic function, and the information coming from other agents. The dynamical evolution of the nth agent is given byẋ where φ n (·) is a nonlinear function, x n corresponds to the vector of state variables, w the vector of environmental inputs, u n the vector of management inputs (i.e., irrigation), and Θ n is the set of parameters related to the nth agent, including soil, crop (e.g., type of cultivar), and management parameters. These parameters can be either constant or time varying. In general, they are known according to a probability distribution, e.g., the soil parameters are stochastic variables accounting for the soil variability. In this study, the dynamical model underlying the crop-soil agents is built up as a combination of elements of three models. A water mass balance proposed by Ritchie [19,20] (CERES-Wheat), the drought index equation proposed by Woli et al. [21] (ARID model), a biomass growth model presented by Zhao et al. [22] (SIMPLE model), and an original waterlogging stress proposed by the authors (the excess of water reduces the crop transpiration and therefore, the biomass production). This set of equations allows computing the biomass produced by every agent as a response to the environmental inputs and the water content in the corresponding patch of soil. The model is developed in discrete-time using a sampling time of 1 day. This is motivated by two reasons: (i) the growth process at the crop scale is assumed to be affected only by the time course of temperature, which enables the use of thermal time (with units degree days, which represents the cumulative time integral of temperature) to compute the daily aging of the crop; and (ii) the environmental inputs are available daily. Then, the discrete-time equations for every agent n = 1, ..., N, can be summarized as follows: x n 2 (k + 1) = x n 2 (k) + h n 1 (k) x n 3 (k + 1) = x n 3 (k) + θ n 11 (1 − h n 2 (k)) + θ n 12 (1 − h n 3 (k)) (4) x n 4 (k + 1) = x n 4 (k) + θ n 13 h n 6 (k)h n 7 (k)h n 8 (k)g n (k)w 4 (k), where the state variables are the water content of every land patch (x n 1 ), the cumulative temperature (x n 2 ), the cumulative temperature until maturity to reach 50% radiation interception due to leaf senescence (x n 3 ), and the biomass (x n 4 ). The fluxes related to the water balance in soil are the crop transpiration ( f n 1 ), the surface runoff ( f n 2 ), the deep drainage ( f n 3 ), and the flux coming from a neighbour ( f n 4 ). The independent inputs include the precipitations (w 1 ) and the solar radiation (w 4 ). The management input is the amount of water irrigated in every patch (u n ). The functions h n (k) are continuous or piecewise continuous functions of environmental inputs (w j ) used to compute the stress factors, and g n (k) is the growth function. Parameters θ n 11 , θ n 12 , and θ n 13 represent the maximum daily reduction in x n 3 due to heat stress, the maximum daily reduction in x n 3 due to drought stress, and the radiation use efficiency, respectively. A complete list of variables and parameters is given in Tables A1 and A2, and the representation of an agent is shown in Figure 3.  Table A1.
Hereafter, the superscript n of all internal functions ( f n (k) and h n (k)), state variables, and parameters is omitted to simplify the notation, keeping in mind that each equation is specific to each agent. The flux components f n (k) related to Equation (2) are computed following the sequence proposed by [23]. The crop transpiration f 1 (k) is given by where θ 1 is the water uptake coefficient, θ 2 is the wilting point, θ 5 is the root-zone depth, and w 2 is the reference evapotranspiration. The surface runoff f 2 (k) is computed as where θ 3 is the initial abstraction. The deep drainage f 3 (k) is estimated as where θ 4 is the drainage coefficient, and θ 6 is the field capacity. The incoming flux f 4 (k) is computed as the sum of outflows f out (k) from all the neighbours, which have a normalized elevation γ higher than the considered agent. The outflow of such an agent is given by where N r is the number of receiving neighbours. The function h 1 (k) in Equation (3) represents the daily mean temperature added to the state variable x 2 , i.e., where w 3 (k) is the mean air temperature and θ 7 is the base temperature for phenology development and growth. The cumulative temperature required to reach 50% of radiation interception during canopy senescence x 3 is increased by heat stress (i.e., faster canopy senescence) and drought stress. The heat stress h 2 (k) is computed as follows where θ 9 is the threshold temperature to start accelerating senescence from heat stress, θ 10 is the extreme temperature threshold when radiation use efficiency (RUE) becomes 0 due to heat stress, and w 6 (k) is the maximum daily temperature. The drought stress h 3 (k) is built based on the ARID index [21], where θ 14 is the sensitivity factor to RUE, and One of the novel features of the model is to consider the negative effect of water excess. Hence, a waterlogging stress factor h 5 (k) is proposed in Equation (14), i.e., where θ 15 is the stress time lower limit, and θ 16 is the stress time upper limit. These limits are specific for every crop, and the values can be retrieved from literature. This index is built based on the anaerobic root stress reported by [24,25]. The rationale behind this index is that the negative effect is proportional to the amount of days (S t (k)) under anaerobic stress. The limits are taken from the literature. The stress shape is also suggested by Steduto et al. [8]. The counting of days under waterlogging stress is computed as The combined negative effect of water excess, lack, and excessive heat is summarized into h 6 (k) using a minimum function due to the rationality that the crop never suffers simultaneously of more than one stress related to transpiration, in other words, The impact of low temperatures on biomass growth rate is calculated as follows where θ 8 is the optimal temperature for biomass growth. Then, the CO 2 impact in biomass growth rate h 8 (k) is estimated as where θ 17 is the relative increase in RUE due to a higher CO 2 concentration, and w 7 (k) is the atmospheric CO 2 concentration. Finally, the growth function g n (k) is estimated as where θ 18 is the cumulative temperature requirement from sowing to maturity, θ 19 is the maximum fraction of radiation interception that a crop can reach, and θ 20 is the cumulative temperature requirement for leaf area development to intercept 50% of radiation.

Irrigation Agent
This agent is the only one that interacts with all crop-soil agents. Its mission is to distribute water for irrigation, either following an arbitrary policy determined by some heuristics or an optimal policy resulting from the formulation of a performance index and the solution of an optimal control problem. There is an array of optimal control policies that could be considered including crop yield, water consumption, and constraints on water availability or distribution. Here, a very simple problem is considered at first, where the performance index corresponds to the average daily biomass, subject to the constraint of maximum available water on a daily basis W max , i.e., where N is the number of agents,x n 4 (k) is the biomass estimate and u n (k) is the irrigation for the crop-soil agent n at day k, respectively.
After the definition of the model agents, it is necessary to examine the flow of information from the sensor.

Sensor
The sensing technology considered in this work is based on images. Images have the advantage of capturing the whole cropping area and can be calibrated with punctual measurements of climatic variables either by a meteorological station in situ, or by ground sensors. For every agent n = 1, ..., N, only x n 1 (k), x n 2 (k), and x n 4 (k) are measurable accordingly to [26][27][28]. The sensor equation is where z n (k) is the vector of sensor outputs, H n is the measurement matrix, and v n (k) is a vector that represents the measurements errors. In this preliminary study, the measurement noise is assumed to follow a normal distribution with covariance matrix R n (k), but other distributions could be considered as well. All measurements are assumed uncorrelated. For the case of x n 4 (k), the measurement is a fraction of the crop biomass (i.e., above ground biomass) [5]. With these signals, the next step is to infer the unmeasurable state variables using state estimation techniques.

Estimator
There are two reasons to incorporate an estimator in the model. First, not all state variables are measurable, and second, it is necessary to incorporate sensor data to adjust the model predictive capabilities. In view of the assumed Gaussian distribution of the noise, a classical solution is provided by an extended Kalman filter (EKF) in discrete-time form [29]. In the framework of the ABM model, the filter is deployed in a decentralized way by associating to each agent a simplified version of the filter where the coupling between agents is partly neglected.
The discrete-time state space representation of a crop-soil agent (Equation (1)) becomes with the consideration of a zero-mean Gaussian process noise ω n (k) with covariance matrix Q n (k). The state space equation is linearised along the estimated trajectory and the Jacobian matrix is obtained by partial differentiation where the coupling between agents is neglected so as to achieve a simple local evaluation The recursive algorithm of the filter begins with the initialisation with the initial statê x n (0), and the error covariance matrix P n (0). Hereafter, the n super index is omitted for simplicity and prior estimates are denoted with an index − . The algorithm is divided into a prediction (or time update) step: and a correction (or measurement update) step using the correction gain K(k): Every agent daily iterates the sequence of Equations (25) to (29). The considerations about the selection and computation of matrices Q(k), R(k) and P(k), follow the recommendations of Reif et al. [30] to guarantee that Φ(k) is bounded and Φ(k) is nonsingular. These recommendations include a positive upper and lower bound for P(k), and a positive lower bound for Q(k) and R(k) (which are positive definite matrix defining the process and measurement noise covariances).
As the exchange flux f 4 is not considered in the Jacobian evaluation, the estimator may be affected particularly regarding the estimation of the water contentx 1 . An ad-hoc procedure is therefore suggested, which consists in "retuning" the variance of the process noise on the first state variable as follows where α is a weight factor and V[w 1 (k)] is the variance of rainfall. There is indeed a direct relation between the rainfall pattern and the intensity of the exchange flux f 4 . This modification allows modulating the uncertainty related to the lack of exchange flux in the model linearisation. An example of tuning of the EKF is given in Appendix B, which is not meant to be optimal, but gives reasonable results.

Agent Spatial Environment
The construction of the ABM requires the selection of the crop and geographic location, e.g., species, cultivar, soil properties. Based on aerial images of the field, image processing is performed, e.g., color to grayscale, binarization, edge detection and region extraction [31,32]. This processing yields a 2-D map, which can be discretized using a spatial grid. The grid size is established based on the properties of the image and the knowledge about the field, e.g., information of the previous crop season. Then, an assignment of the N agents is carried out. The minimum ground area for parameter calibration is equal to 1 m 2 , which is selected as the lower limit of agent size. Then, the maximum number of agents is the ratio between the minimum ground area and the total crop area. Once a label is assigned to an agent, the algorithm extracts extra information from the image, e.g., the normalized elevation parameter γ for every agent [33,34]. Subsequently, a graph is created, as well as the adjacency matrix and the matrix of neighbours. As the number of agents increases, the graph can be used to decompose the cropping region into smaller portions while keeping the relationships between neighbours.

Model Implementation
The survey paper [35] presents a comprehensive comparative literature review of the state-of-the-art software development in agent-based modelling. Mesa, a modular agent-based modelling framework in Python [36], appears as an efficient tool to develop the crop-soil model proposed in this study. It allows users to create agent-based models using built-in core components, such as spatial grids and agent schedulers, or customized implementations, to visualize them using a browser-based interface, and to analyse the results using Python data analysis tools. The modules are grouped into three categories:

1.
Modeling: Modules used to build the models themselves: a model and agent classes, a scheduler to determine the sequence in which the agents act, and space for them to move around.

2.
Analysis: Tools to collect data generated from the model, or to run it multiple times with different parameter values.

3.
Visualization: Classes to create and launch an interactive model visualization, using a server with a JavaScript interface.
Once the image processing routine is executed and data is stored in text files, the model is built. This model is split into three files. The first file contains the model definition, where the irrigation agent is incorporated as a "method" (an object-oriented programming procedure) of the model, as there is no physical location of this agent (see Algorithm A2 of Appendix C). The second file contains the crop-soil agent definition and methods (see Algorithm A1 of Appendix C). Finally, the last file contains the web visualization and the plotting routine for the state variables. The whole model is intended to advance daily (k = 1d). However, every agent can behave asynchronously by setting its own self_time variable.

Case Studies
In order to assess the functionality of the model, three different scenarios with heterogeneous soil are shown in Figure 4. These examples only consider the water supply problem, but no nutrient stress. The sample fields are located in the region of Boyacá For all scenarios, the crop is wheat of the same cultivar (i.e., yecora rojo) and the land size is 1600 m 2 . The study area is a square shaped land, and no incoming and outgoing fluxes beyond this area are considered. After image processing, N = 16 agents are chosen, and the parameters Θ n are adjusted to match the regional reported values. All parameters are listed in Table A2. The first scenario is a flat terrain (i.e., slope less than 3%) in which the root-zone depth θ n 5 changes according to a uniform probability distribution as suggested in Woli et al. [37]. This parameter was selected after evaluating the impact that the parameters θ 1 to θ 6 have on the state variable x 1 , following the methodology proposed by Woli et al. [37]. This study shows that the root-zone depth parameter has a 42% influence on the output based on the Fourier amplitude sensitivity test. This scenario is evaluated during the two sowing periods each year. The second scenario considers a sloping ground (i.e., slope around 20%) with the same climatic data as in the previous scenario. Finally, the third scenario is concerned with a rugged ground that allows evaluating different interactions between the agents due to the topography. The soil is mostly loamy and sandy-loamy in all scenarios.
For this case study, aerial images are only considered to configure each of the scenarios and establish the initial conditions of the parameters and variables. Subsequently, the daily information from the sensors is generated by the model (i.e., Equations (2) to (5)) with the addition of Gaussian noise with zero mean and standard deviation in the range described by [26,27], and [28] for the variables x 1 (k), x 2 (k) and x 4 (k), respectively.

First Scenario
In this scenario, there is no exchange of water between agents. Instead, the root depth zone of every agent θ 5 is chosen from a uniform distribution U(300, 600). The assignment of agents is presented in Figure 5. The value of x n 1 (0) is at field capacity for all agents. To assess the irrigation effect on every agent, the model runs first with no irrigation and then with the irrigation policy. For illustration purposes, only agent 1 is considered. All state variables are displayed in Figure 6 as all of them could be considered to make management decisions, i.e., x 1 for water use, x 2 to control the phenological development of the crop, x 3 to evaluate the effect of stress factors in crop transpiration, and x 4 to evaluate production. The estimator is working satisfactorily with a convergence time of about 30 days. x 3 , which is the only unmeasured variable, is estimated with a slight offset. The error of estimation of the biomass x 4 , which is used in our performance index to determine the irrigation policy, is ranging between 0.5 and 5% in the several simulation tests, see Table 1.  Figure 7 shows the normalized soil water content with respect to the field capacity without and with irrigation. In the latter case, the water level is around field capacity. Moreover, the water stress indexes (i.e., drought and waterlogging) indicate that only if the field capacity is exceeded for more than 3 days the waterlogging stress is triggered. After irrigation, the drought stress decreases and the field becomes prone to waterlogging stress when the precipitations increase. The positive effect of irrigation is seen in the increment of the production of biomass. Table 1, shows that the water quantity used for irrigation is lower than the total rainfall except in the situation of dataset 2 where it is 132 mm higher.

Second Scenario
With a sloping ground, the water flux f n 4 (k) = 0 from the most elevated side of the terrain as represented by the arrows in the graph (see Figure 8). To illustrate this situation, agents 1 and 16 are selected, which are in the highest and lowest part of the ground, respectively, and also are at maximum distance from one another.
The trajectories of Agent 1 and 16 and their state estimates are shown in Figure 9 for state variables x 1 , and x 3 . The performance of the state estimator is satisfactory for agent 1, but the estimation of the water content x 1 of agent 16 suffers from transient errors particularly in periods of strong rainfall, where the water content is underestimated. The tuning parameter α can be used to partly alleviate this situation, but is not enough to guarantee a good estimation in all cases. After irrigation, state x 3 decreases as a result of the reduction in the drought stress and the increase of the crop transpiration.  Figure 10 shows the water effect on biomass production for the first and fourth climatic datasets.
The waterlogging stress during the first 30 days of growth has no significant impact on the final biomass production. However, after this period the irrigation increases the waterlogging stress of agent 16 under strong rainfall, and delays the biomass production. The biomass production x 4 increases by 22, and 29% for agent 1 but decreases by 2.4, and 14% for agent 16 for datasets 1 and 4, respectively (see Table 1).

Third Scenario
In this scenario there are several interactions between agents ( Figure 11). Agent 2 represents the highest point in the field and agent 10 corresponds to the lowest location.
The estimation of the state variables x 1 and x 3 with and without irrigation for the first meteorological dataset is shown in Figure 12. Again, it can be observed, that the water content of agent 10, which is at a low point and therefore receives more water, is underestimated under stronger precipitations. The estimation of x 3 with and without irrigation is subject to a slight offset. The reduction of x 3 by irrigation is related to the decrease in the drought stress and therefore an increase of crop transpiration. Figure 13 shows the water effect on biomass production for datasets 1 and 4. For the first dataset it is noticeable that the gusts of rain generate a waterlogging stress even for the agents in the highest points of the terrain (e.g., agent 2). This effect is enhanced by irrigation. Moreover if the waterlogging stress is triggered in the growth stage, the biomass production is affected negatively. The positive effect of irrigation on the biomass production is about 30% for agent 2 and between 1.2 to 2.5% for agent 10 for datasets 1 and 4 (see Table 1). In order to compare the performance of the irrigation policy in the various cases, a consume ratio is proposed as where W max = 20 mm, N = 16 and M = 160 (the harvest day is t f = 160 d). Table 2 presents a global view of the biomass production and water consumption for irrigation. The latter represents 6 to 11% of the maximum available water in all cases. The positive effect of irrigation is evidenced by a global rise of the production ranging between 10 to 62% depending on the case (land topology and weather). Obviously, the total biomass maximization does not necessarily increase the individual agent biomasses as shown in Table 1, where for scenarios 2 and 3 the agents that are located at low points (i.e., agents 16 and 10) suffer from waterlogging instead of drought. The estimator is performing reasonably well, with usual errors in a range of 5%, but with transient discrepancies up to 20% for low points of uneven lands and strong precipitations.

Conclusions and Perspectives
An agent-based model of an agricultural system is presented, which can be used to assess different irrigation schemes and production scenarios. The dynamic crop model underlying each agent is built based on existing models, but includes a novel waterlogging stress factor, proposed to consider the negative effect of water excess in soil. The ABM also includes an estimator of the unmeasured variables in the form of an extended Kalman filter, and a simple optimal irrigation strategy based on the average daily biomass. This framework is flexible and scalable, which allows the user to add new agents (or cropping areas) with diverse crops and different soil properties. Furthermore, crop-soil agents can be resized to catch the main soil changes across the cropping land based on soil properties and topographic variations. The ABM model is coded in the Python environment Mesa which is suitable for this purpose. However, there is the need to develop complementary modules to enhance the model functionality, e.g., GIS integration.
Future research entails the continued validation of the proposed model, for instance following the methodology of [38][39][40][41], including the use of global sensitivity analysis to assess the influence of the key parameters and possibly simplify the model structure in order to ease parameter estimation. Parameter identifiability is a critical aspect to develop a reliable prediction model, based on the collection of experimental data.
In addition, the consideration of other estimation schemes should be the subject of further investigation. Indeed, the current study was limited to a simple decentralized extended Kalman filter scheme, which shows poor performance for uneven grounds where water can accumulate at low points. Future developments therefore entail the use of centralized estimation schemes (i.e., encompassing all the agents), and the consideration of alternative estimation methods in order to avoid model linearisation and to handle non-Gaussian noise distributions as well as measurement correlations, which might be more appropriate in the extraction of information from image analysis.
Finally, the definition of the irrigation agent might be refined to describe the situation where irrigation is not available everywhere and/or not in the same supply. Different performance indexes might also be investigated, including weather forecast, in the framework of model predictive control. After validation and additional tests, an updated version of the ODD protocol [42] will be made available to the readers upon request.

Data Availability Statement:
The data presented in this study is available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.  Threshold temperature to start accelerating senescence • C 34 [22] from heat stress θ 10 The extreme temperature threshold when RUE becomes 0 • C 45 [22] due to heat stress θ 11 The maximum daily reduction in x 3 due to heat stress • C d 100 [22] θ 12 The maximum daily reduction in x 3 due to drought stress • C d 25 [22] θ 13 Radiation use efficiency (above ground only g MJ −1 1.24 [22] and without respiration) θ 14 Sensitivity of RUE to drought stress (ARID index).

Appendix B. Jacobian and EKF Tuning
The Jacobian is computed as follows where Heat stress by Equation (11) # High temperature 25 Temperature stress by Equation (17) # Low temperature 26 Drought stress by Equation (12) from Equation (13) # Lack of water 27 Waterlogging stress by Equation (14) from Equation (15)