A disaggregated system dynamics and agent-based modeling of the water-energy-food nexus for optimizing water allocation

This study introduces an integrated modeling framework rooted in the WEF nexus to assess and optimize agricultural resource allocation in the Zayandehrud River Basin. The methodology consists of three primary components: (1) simulation of the WEF system using a hybrid SD-AB approach, (2) multi-objective optimization using NSGA-II, and (3) prioritization of Pareto-optimal strategies through the AHP as shown in Fig. 1.

Flowchart of the simulation–optimization model.

SD-AB simulation model

The SD–AB hybrid model captures both the macroscopic feedback across water, energy, food, and economic subsystems and the behavioral diversity of decentralized agents (e.g., farmers, water users, subregions). SD modeling is applied to represent dynamic feedback and flows between WEF components over time, incorporating causal loop diagrams to structure reinforcing and balancing mechanisms^7,18.

Each of the 14 agricultural subregions within the AB part of the SD–AB modeling framework is characterized as an autonomous agent. The agents are heterogeneous in terms of groundwater availability, crop type, irrigation efficiency, and economic conditions. The behavior of each agent—expressed by their water use, field area decisions, and the associated costs and benefits—is modeled separately for each subregion. The results produced by each of the agents are synthesized within the system dynamics model in order to determine their collective effect on basin-scale variables, including the surface water flux and the groundwater balance. There are no direct interactions among the agents, but their activities are coupled via shared resource dynamics and system-level feedback mechanisms.

The schematic representation of the interactions among the internal components of the WEF nexus system is shown in Fig. 2. This diagram is the basis for SD modeling, serving a crucial role in visualizing and capturing the dynamic behavior of the WEF nexus. The functionality of the WEF nexus is thoroughly investigated using the SD model, which allows for an in-depth exploration of the interrelated factors within the system.

The Zayandehrud River basin. This map was created by the authors using ArcMap version 10.7.1 (Esri, https://www.esri.com/en-us/arcgis/products/arcgis-desktop/overview)”.

The CLD shown in Fig. 3 delineates the behavior of the modeled WEF system, identifies potential feedback mechanisms, and provides insights into the long-term impacts of various policies and interventions. The model subsystems are interconnected through nexus variables, which represent the interactions and interdependencies among the water, energy, food and economics subsystems. These nexus variables are essential for analyzing how changes in one subsystem can influence the others and for understanding the overall system dynamics¹⁹.

Schematic of the interactions that exist in the water-energy-food-economic system elements (SW: surface water, GW: groundwater).

Water availability -shaped by agricultural water demand and groundwater abstraction- is identified as a nexus variable impacting the food subsystem, particularly through irrigation practices, and the energy subsystem, notably through hydropower generation. Fluctuations in water availability include cascading effects throughout the energy and food subsystems, creating significant feedbacks and interconnections²⁰.

The food subsystem is intricately linked to the water subsystem through multiple key pathways, as illustrated in Fig. 3. It relies on the water subsystem to fulfill crop water requirements. The total water supply is connected to agricultural water demand, which subsequently affects crop yield. As crop yields increase, the demand for irrigation water may rise correspondingly, highlighting the dependence of food production on adequate water resources²¹. Furthermore, the relationship between population growth and food demand illustrates how increasing food requirements exacerbate reliance on the system’s water resources. In this context, “Cultivated land” emerges as another significant nexus variable that connects the water and food subsystems, reflecting relationships among these components.

The water subsystem

The water subsystem includes surface water (SW), groundwater (GW), water demands (WD), and return water to water resources (InfSWD and InfGWD).

The surface water balance is presented in Eq. (1). All dimensions are expressed in million cubic meters (MCM) and time is expressed in years.

in which, \({SW}_{net}\), SW, P, Q, TW, E, InfSW, and \({SWD}_{Total}\) denote respectively the surface water balance, available surface water, rainfall, inflow (runoff), transferred water, evaporation, surface water infiltration into groundwater (i.e., recharge), and total surface water demand at time t. SWInfCo is the coefficient of surface water infiltration into groundwater.

The groundwater balance is presented in Eq. (3) and infiltration from rainfall is calculated using Eq. (4).

in which \({GW}_{net}\), GW, InfP, InfSWD, InfGWD, and \({GWD}_{Total}\) represent the groundwater balance, available groundwater, infiltration from rainfall, return flow from surface water usage, the return flow from groundwater usage, and the total groundwater demand at time t, respectively. PrecipInfCo and Precipitation_t denote, the precipitation infiltration coefficient and the amount of precipitation, respectively.

InfSW denotes surface water infiltration into groundwater, as presented in Eq. (5).

in which SWD_Agr, SWD_Ind, and SWD_Dom, denote agricultural, industrial, and domestic surface water demand, respectively, and InfCoSWD_Agr, InfCoSWD_Ind, and InfCoSWD_Dom denote the coefficients for agricultural, industrial, and domestic surface water demand, respectively.

Return flow from groundwater usage is calculated as Eq. (6).

in which GWD_AgrC, GWD_AgrH,, GWD_Ind, and GWD_Dom denote, respectively, the groundwater demand for agricultural crops, agricultural horticulture, industrial use, and domestic use, and InfCoGWD_AgrC, InfCoGWD_Ind, and InfCoGWD_Dom denote, respectively, the coefficients for groundwater demand for agricultural crops, agricultural horticulture, industrial use, and residential use.

The infiltration values from runoff in plains and heights were obtained using Eqs. (7) and (8).

in which RunoffPlainCo and RunoffHeightCo denote, respectively, the conversion coefficients of precipitation to runoff in plains and heights (elevated areas), and InfCoRunoffPlain and InfCoRunoffHeightCo denote, respectively, the runoff infiltration coefficients.

Initial parameter values were determined based on historical hydrological records (2001–2021), supplemented by values reported in the literature and adjusted through calibration using observed surface water and groundwater levels. Table 1 lists the model parameters of fundamental significance, both their preliminary estimates and final values, thereby illustrating the adjustments achieved through the process of model calibration.

The sum of the inputs to and outputs from the aquifer expresses the changes in the storage volume [Eq. (9)].

in which \(\Delta\) V, S, A, and ΔH represent the storage change in the aquifer, the storage coefficient of the aquifer, the area of aquifer (square kilometers), the fluctuation of groundwater level (meter), respectively. If there is a balance between inputs and outputs the storage change would equal zero.

The storage coefficient is a hydrodynamic characteristic, which varies with the specific storage factor and the thickness of the aquifer, varies between 5 × 10^–5 and 5 × 10^–3. The storage coefficient of aquifers in the areas is listed in Table 1.

The surface water demands and groundwater demands are presented in Eq. (10) and Eq. (11), respectively

in which \({SWD}_{Agr}\) \({SWD}_{Ind}\), \({SWD}_{Dom}\), \({SWD}_{Env}\), denote respectively surface water demand ((MCM) in the agricultural, industrial, domestic, and environmental sector, and \({GWD}_{Agr}\), \({GWD}_{Dom}\), \({GWD}_{Ind}\), denote respectively groundwater demands (MCM) in the agricultural, industrial, and domestic sectors.

Domestic water demand and total agricultural water demand is calculated using Eqs. (12) and (13), respectively.

in which \({Q}_{Dom}\) and N denote respectively domestic water consumption per person and population size.

in which \(m\) is an index denoting a crop type (there are \(M\) crops considered as decision variables in this work).

The Food Subsystem

The food subsystem estimates crop water requirements using FAO Penman–Monteith²² evapotranspiration models, integrating cultivation pattern, yield, irrigation efficiency, and cropping patterns (Eqs. (14) to (18)).

in which I_n, ET_o, P_e, G_e, K_c, I_g and E denote, respectively, the crop net water requirement, actual crop evapotranspiration, effective rain, groundwater absorbed by crops, crop coefficient, gross water requirement, irrigation efficiency. R_n, G, e_a, e_d, T,\(\gamma\), Δ, U₂ denote, respectively, the net radiation flux (MJ/m²day), soil heat flux (MJ/m²day), mean saturation vapor pressure (kPa), mean ambient vapor pressure (kPa), air temperature at 2 m height, psychrometric constant (kPa/°C), slope of the saturation vapor pressure vs air temperature function at the air temperature at a height of 2 m, the wind speed at a height of 2 m (m s⁻¹).

The energy subsystem

The energy subsystem calculates energy consumption primarily for groundwater extraction based on hydraulic head (\({lift}_{t}\)) and water mass (\({mass}_{t}\)), establishing water–energy coupling^23,24 in Eq. (19).

In which, \({E}_{t}\) and \(\varepsilon\) denote, respectively, the energy requirement in period t and the efficiency of the energy generating system (typically between 0.7 and 0.9).

The economic subsystem

The economic subsystem computes agricultural net benefits using crop-specific costs, prices, and yields, connecting profitability to WEF performance indicators through Eqs. (20) and (21), respectively:

in which \({CY}_{j}^{t}\), \({Production}_{j}^{t}\), \({Area}_{j}^{t}\), \({NetBen}_{j}^{t}\), \({Price}_{j}^{t}\), and \({Cost}_{j}^{t}\) respectively denote the yield of crop \(j\) in period t (Kg ha⁻¹), production of crop \(j\) in period t (Kg), cultivation area of crop \(j\) in period t (ha), net benefit from crop \(j\) in period t (Rial), price and cost of crop \(j\) in period t (Rial).

The AB component models the agricultural subregions as autonomous agents with the AnyLogic simulation software²⁵. Each agent responds dynamically to simulated groundwater levels, profitability, and water access, influencing spatial allocation decisions. A detailed description of decision variables is found in Sect. “Multi-objective Optimization”.

Data sources

This study’s WEF simulation uses 20 years (2001–2021) based on the datasets listed below, with policy relevance and historical calibration listed in Table 2. The model operates with a yearly time step, which is consistent with the temporal resolution of available meteorological, hydrological, and agricultural, energy, and socio-economic data and enables the analysis of long-term trends in water, energy, and food system dynamics.

The model uses the following datasets:

• Hydrological data: Precipitation, evapotranspiration, infiltration, surface/groundwater flows and returns;

• Agricultural data: Crop-specific land use, yield, irrigation need;

• Energy data: Electricity use for irrigation, with 85% allocated to agriculture;

• Economic data: Crop prices, production costs, per capita food demand, population.

The last column of Table 2 offers a comprehensive view of the principal sources of primary data and their corresponding agencies in charge of each kind of dataset. The agencies are the Iran Water Resource Management Company (IWRMC), Agricultural Jihad Organization of Isfahan Province, Ministry of Energy of Iran, Ministry of Agriculture (Jihad), and the Statistical Center of Iran.

Missing data in hydrological and meteorological data sets were filled in by linear regression using data from nearby stations. Missing data in agricultural and energy-related data sets were completed with national reports and planning documents. These actions assured completeness and consistency of input data for the 2001–2021 modeling horizon.

Model calibration and validation

The water subsystem calibration was performed using 70% of the historical data (13 years), while 30% (6 years) was used for validation.

A combination of graphical techniques, dimensionless statistics, and error indices was employed to evaluate model performance²⁸. Specifically, the correlation coefficient (R), Normalized Relative Root Mean Square Error (NRMSE) and, Nash–Sutcliffe efficiency (NSE). Performance ratings for the NRMSE and the NSE values are listed in Table 3.

The SD-AB model demonstrated excellent performance with respect to surface water simulation with a calibration NSE of 0.988, R = 0.996, and NRMSE = 0.034. During the validation period, NSE remained high at 0.933, with R = 0.966 and NRMSE = 0.106, as listed in Table 4. The SD-AB model accuracy with respect to groundwater simulation varied across subregions and generally indicated reliable performance. For instance, in sub-basin 4205, the validation NSE reached 0.91 with R = 0.98, while sub-basin 4201 exhibited lower performance during validation (NSE = 0.67, R = 0.93). In the subregions where there are the lower calculated values for the NSE in the validation period the lower performance is due to input data of poorer quality or to local hydrologic heterogeneity not completely represented by the aggregated model structure. In addition, the accuracy of the model is diminished in situations of uncontrolled human intervention—such as over-abstraction of groundwater—or discrepancies in the observational records.

The full range of calibration and validation results across subregions is listed in Table 5. The model’s predictions aligned closely with the observed trends and values, demonstrating a strong agreement between the simulated and actual data. It is important to note that groundwater levels for regions 4210, 4215, and 4216 are not recorded.

Multi-objective optimization

Sustainable strategies under competing objectives were identified by coupling the model with the NSGA-II evolutionary algorithm²⁹. The optimization sought to:

• Minimize water consumption (Z₁) (defined by Eq. (13)),

• Minimize energy consumption (Z₂) (defined by Eq. (19),

• Maximize net agricultural benefit (Z₃) (defined by Eq. (21)).

The model incorporates nine decision variables, which are the agricultural surface water uses (X₁), and groundwater allocation for wheat (X₂), barley (X₃), maize (X₄), alfalfa (X₅), rice (X₆), potato (X₈), other crops (X₉), and orchards (X₉).

The NSGA-II algorithm is configured with a population size of 20, crossover rate of 0.2, and 50 generations, balancing computational efficiency with solution diversity³⁰. These values were determined through trial-and-error testing, balancing the model’s computational demands with the need for adequate solution diversity. Larger population sizes or more generations only slightly improved the quality of the solutions while greatly increasing the computational load.

Several constraints were imposed in the optimization model. Limitation of surface water (Eq. (22)) and groundwater (Eq. (23)): The water allocated for different uses cannot exceed the available water resources (SW and GW).

Limitation of the water balance (Eq. (24)): Establishing the balance equation is one of the conditions that must be met in the study region, and it is calculated based on evapotranspiration and precipitation.

Limitation of the area under cultivation (Eq. (25)): The cultivated land (\({Area}_{t}^{j}\)) for crop j in period t (ha) cannot exceed the total arable land (\({Area}_{Total}^{t}\)).

Limitation of food security³¹ (Eq. (26)): Food security means the access of all people at any time to sufficient, healthy and nutritious food that is necessary to meet nutritional needs and lead an active and healthy life³². The food demand (FD) for each person is the same as the per capita consumption of each person per year, which is determined by the total of drink production and imports per population.

The WEF nexus multi-objective optimization model was solved with the NSGA-II optimization algorithm^{29,33,34,35,36,37,38,39}.

solution evaluation viaAHP

The AHP method is used to rank alternatives based on criteria defined by stakeholders given the trade-offs among Pareto solutions: water availability, environmental impact, and cost-effectiveness^40,41. The most appropriate strategy is selected through pairwise comparisons and weighted scoring and aligning technical performance with policy priorities.

Case study

The Zayandehrud Basin (Fig. 4) in central Iran covers ~ 27,000 km² and is an important critical region for agriculture, industry, and urban supply. Persistent water conflicts, high agricultural demand, and ecological degradation in the Gavkhuni wetland make this basin an ideal setting for WEF-integrated modeling^42,43.

The developed conceptual loop diagram (CLD) model’s graphical representation.