The simulated environment consists of a 210 × 210-pixel field representing a 210 × 210 m section of the larger experimental plot. The spatial arrangement mirrors the real-life design of the integrated agroforestry system, consisting of I. paraguariensis (“yerba mate”) plants as the base crop and nine consociated tree species. I. paraguariensis plants are positioned on a 20 × 20 grid, with coordinates at {10, 20, 30, …, 200} in both x and y dimensions. Tree agents are placed on a parallel grid offset by 5 pixels, with x-coordinates {5, 15, 25, …, 205} and y-coordinates {5, 15, 25, …, 205}. Each of the nine tree species occupies two adjacent columns, ensuring that every I. paraguariensis plant is consociated with four equidistant tree neighbors. Conversely, each tree interacts with four surrounding I. paraguariensis plants. A control area, devoid of trees but containing I. paraguariensis, occupies the final two columns (x = 195 and x = 205).

In summary, the model initially contains 400 I. paraguariensis plants arranged in a 20 × 20 grid. Of these, 360 plants are consociated with trees, while 40 plants, located in the control area (two rightmost columns), grow without tree association. The simulation includes 42 trees per species, for a total of 378 tree agents representing nine different species. Collectively, these agents occupy a 200 × 200 m section of the larger 10-hectare experimental field, which contains approximately 7,460 trees and 37,000 I. paraguariensis plants in total (Comolli et al., 2024, 2025a). The tree species included in the model are listed in Table 1, and representative visualizations of the simulated field under different initial soil conditions are shown in Supplementary Figure 1. The spatial configuration mirrors that of the real-world experimental design, ensuring that model-derived dynamics correspond directly to measurable field processes.

The design of the real-life experimental plot, the selection of tree species, and the allometric parameters used for calibration are detailed in (Comolli et al., 2024, 2025a). For clarity within the model’s code, specific abbreviations are used for each species: I. paraguariensis (ilex); Peltophorum dubium (Cana); Handroanthus heptaphyllus (La); Cordia trichotoma (Pet); Parapiptadenia rigida (An); Araucaria angustifolia (Ar); Balfourodendron riedelianum (Gu); Grevillea robusta (Gr); Paulownia sp. (Ki); Toona ciliata (Toona); Control (Con). The details of the original randomized trial plot, which spans 10 hectares and is divided into four equivalent blocks with varying species distributions, are available in Supplementary Table 1 and in (Comolli et al., 2024, 2025a).

The attributes of the agents in our model are parameterized based on measured field data, as presented in Comolli et al (Comolli et al., 2024, 2025a). Diameters at breast height (DBH —at 1.30 meters from the ground) and heights (Hts) were measured in the field across all tree species and provide estimates for their growth rates presented in Table 2. These data, which are the most commonly used measures of tree growth (Sterck and Bongers, 1998; Sumida et al., 2013), were measured at the ages of 4, 5, and 8 years, and are presented in Supplementary Table 2, taken (with permission) from (Comolli et al., 2024) Species-specific growth rates were derived directly from measured increments in DBH and height, as the average of the two (see “Overview of the NetLogo Code” in Supplementary Materials). This simplification captures the development of the trees during the first few years but cannot be used to model the long-term development of trees. In this model, we focus on a five-to-ten-year period, primarily examining the interplay of competition and synergies between the trees and the harvest-producing plant, I. paraguariensis, as the plants develop from their juvenile growth stages to productive maturity. The measured photosynthetically available radiation (PAR) for each species underpins the modelling of competition caused by shade and physical occlusion incorporated into the model (see “Overview of the NetLogo Code” in Supplementary Materials). Additionally, the model includes human interventions, such as harvesting I. paraguariensis and managing the trees through crown trimming and thinning (reducing the number of trees). The parametrization therefore is based on data acquired over a period of ten years, as published in Comolli et al (Comolli et al., 2024, 2025a).

The simulation commences with the choice of one of five distinct soil quality levels, ranging from degraded to restored, which are represented by a dimensionless ‘energy’ proxy that plants consume for growth. This framework reflects our system’s objective of sustainable soil improvement over time. Each species’ initial energy value is contingent on the designated soil quality level, with the variable ‘et’ (energy trees) assigning dimensionless energy values to all plants based on a soil-factor multiplier. This proxy is a simplified yet realistic representation of soil fertility, calibrated to qualitatively mimic observed field conditions.

As I. paraguariensis plants and trees grow, they draw energy from the soil, thus depleting its quality. Their energy is periodically reduced through harvesting, trimming, and thinning interventions. This dimensionless cyclical process of growth and intervention facilitates the reintroduction of carbon-rich matter to the soil, which the model tracks as dimensionless ‘soil-multipliers.’ These multipliers subsequently update the ‘soil-energy’ and ‘soil’ values (dimensionless). SOM, a critical component acting as both a buffer and a carbon source, also gradually enhances the efficacy of external inputs, particularly synthetic fertilizers. As soil quality changes, plant energy allocations are recalibrated, simplifying the cycle representation while maintaining realism (Bashir et al., 2021; Hoffland et al., 2020; Tian et al., 2023; Abdu et al., 2023; Gao et al., 2011).

The growth rate for all plant and tree agents is species-specific and dependent on soil quality. These growth rates are modelled from allometric field data collected in 2014 and 2015, as presented in Supplementary Table 2 (Comolli et al., 2024). Photosynthetically active radiation (PAR), measured in 2018 (Comolli et al., 2024), serves as a control to verify that simulated competition attributable to consociations scales appropriately with field experimental observations. Roots are not included in this model.

Agent interactions are categorized as either direct or indirect. Direct interactions, such as physical competition and PAR occlusion, directly impact I. paraguariensis plants through shading, reducing their growth rate when light occlusion exceeds 30%. Concurrently, shade promotes soil humidity conservation, which is beneficial. Indirect interactions primarily involve soil enrichment from decomposing trimmed branches, harvest remnants, and thinned logs. These represent endogenous sources of organic matter, which exhibit minor species-specific differences. SOM, in turn, enhances the plant availability of exogenous fertilizers, constituting a second-order, indirect interaction within the model. Additionally, each tree species uniquely affects soil micronutrient composition, such as phosphorus content, further contributing to indirect agent interactions. While challenging to parameterize directly, these complex factors could be incorporated through a maximum entropy formulation.

Human interventions, including I. paraguariensis harvesting and tree trimming and thinning, are integrated as periodic reductions in plant biomass that directly contribute organic matter back to the soil, thus closing the growth-reincorporation cycle. The management of the agroforestry system including the schedule of I. paraguariensis harvest, trimming and thinning of the different tree species as determined by their growth, external fertilizer additions, as well as shade competition, incorporation of SOM from the biomass harvested from the trees and left on the ground are described in the section “Overview of the NetLogo Code” of the Supplementary Materials. The algorithm is included in Appendix I of the Supplementary Materials.

Soil energy is represented as a dimensionless index integrating fertility, organic matter availability, and nutrient-use efficiency. Absolute field measurements were rescaled to preserve relative differences and empirically observed thresholds rather than physical units. This approach maintains internal commensurability while enabling comparison across management scenarios and initial soil states.

The ‘State-of-soil’ plot depicts the temporal evolution of soil fertility, represented as ‘energy’ within the model. The initial soil state is configured using the ‘Fertility’ slider, which categorizes the soil into five distinct types ranging from zero to ten thousand dimensionless units. The five selected soil quality levels, along with the represented energy flow between soil and plants, are dimensionless and qualitative. The critical determinants are that their relative scaling and evolution accurately reflect observed reality, and that this simplified representation facilitates the discernment of system behavior beyond readily apparent observations.

As the simulated system evolves, the soil classification may dynamically change. Tree agents deplete soil energy through nutrient absorption; however, pruning and thinning interventions replenish the soil with organic matter. This process incorporates atmospheric carbon and nitrogen from rainfall (released during electrical storms). These endogenous factors not only enrich the soil but also enhance the efficacy of exogenously added fertilizers. SOM, a critical component serving as both a buffer and a carbon source, gradually increases the solubility and, consequently, the availability of exogenous industrial fertilizers (Comolli et al., 2024, 2025a).

Soil quality is discretized into five empirically grounded ranges based on long-term field observations and published measurements of soil organic matter dynamics in the experimental agroforestry system:

These simulations represent a system established on fully restored soil (‘Soil 1’), typically achieved through intensive restoration efforts such as establishing artificial forests with fast-growing species like eucalyptus or conifers. These conditions match those of the implementation and experimental measurements of the multispecies agroforestry trial, spanning ten hectares, described in Comolli et al (Comolli et al., 2024, 2025a). While this higher initial soil quality enables examination of system trajectories attending optimization from favorable initial conditions rather than restoration, our simulations reveal that even these favorable conditions require careful management. Without any reposition, the system cannot maintain itself, as extraction and plant growth exceed the endogenous fertility maintenance capacity. We examine two management scenarios: fixed minimum repositions and adaptive reposition.

Scenario. We simulate the system starting from a fully restored starting condition (soil quality index = 8,100; “Soil 1”) with a fixed external reposition set to 30% of the extracted biomass and a 40% harvest rate. Despite high initial fertility supporting robust early growth, soil quality first declines for several years until endogenous SOM accrues enough to flip the system into a positive feedback loop. Figure 1 shows representative spatial snapshots across the 10-year run; Figure 2 tracks soil quality and cumulative harvest over time.

Initial phase (0–3 years). Soil quality drops from “Soil 1” to the “Soil 3” range (≈5,000 units) by 137 weeks (~2.6 years) and remains near that level through year 4 (Figure 2a). Plant growth is still robust during this window (Figure 1a), buoyed by the high starting fertility. Fast-growing trees undergo 50% thinning, but incorporation of biomass into soil lags; early soil gains are driven mainly by external inputs. As trimming and thinning begin to contribute measurable SOM, the fertilizer-to-harvest ratio declines, and net soil gains stabilize (Figures 3, 4).

Transition phase (3–5 years). By 4 years 5 months (228 weeks), the fastest-growing trees have already undergone 50% thinning and are regenerating (Figure 1b). The slower species will undergo 50% thinning during this period, along the cutting of the remaining 50% of the fastest growing trees. The largest I. paraguariensis plants occur between recently thinned tree lines—zones of minimal shading and high SOM input from felled biomass. Conversely, the smallest plants appear in the control rows without trees, where soil improvement is limited. A critical transition occurs as accumulated endogenous SOM and stable external inputs jointly trigger a self-reinforcing positive feedback loop, driving soil fertility above initial values (Figures 1b, 2a, 3, 4). The I. paraguariensis harvest, which previously declined during soil depletion, begins a sustained recovery (Figure 2b).

Mature phase (5–10 years). After year five, the system develops increasing structural heterogeneity. Multiple thinning cycles across species generate alternating patches of mature and regenerating trees (Figures 1a–d). The largest I. paraguariensis individuals consistently occur near fast-growing species, benefiting from localized SOM enrichment following thinning. The fertilizer-to-harvest ratio declines steadily as soil quality improves (Figure 4), indicating that fertilizer costs represent a diminishing share of harvest value. By the end of the 10-year simulation (520 ticks), soil quality has increased by ~6,000 units—from ~8,000 to ~14,000—and cumulative I. paraguariensis harvest reaches ~4,500 units. Across all individuals, sizes exceed those observed in simulations starting from lower soil indices, reflecting a self-sustaining cycle of SOM-driven soil improvement, enhanced growth, and reduced external input dependence.

The combined patterns observed in Figures 3 and 4 reveal how the system transitions from dependence on external fertilizers to a self-reinforcing equilibrium driven by endogenous SOM accumulation. Initially, external inputs dominate the nutrient balance, sustaining productivity while soil reserves are rebuilding. As thinning, trimming, and harvest residues progressively enrich the soil, endogenous SOM becomes the primary driver of fertility gains, leading to a steep rise in soil quality (Figures 3b, 4a). This shift coincides with a marked increase in total biomass production and a sharp decline in the fertilizer-to-harvest ratio (Figure 4b). By year ten, the system operates in a low-input, high-yield regime where most nutrient cycling is internally maintained. The efficiency gain—expressed as declining fertilizer costs relative to total harvest—demonstrates the establishment of a virtuous cycle linking soil restoration, biomass accumulation, and long-term productivity.

This scenario examines an adaptive repositioning strategy that begins at 25% of extraction and automatically adjusts inputs in response to harvest feedback. Simulations start from fully restored soil (soil quality index = 8100, Soil 1) with a 40% harvest rate. Similar to the fixed‐reposition case, the system initially exhibits a transient decline in fertility before endogenous SOM accumulation and feedback-regulated inputs jointly drive recovery. This gradual buildup initiates a positive feedback loop that elevates soil quality beyond its initial state and sustains higher annual harvest volumes (Figures 5–8). Although the high initial fertility temporarily declines—producing reduced early harvests (Figures 6a, b)—adaptive inputs, continuously calibrated to match extraction (Figure 7), progressively enhance biomass production (Figure 8a). By approximately year five, the cumulative effect establishes a self-reinforcing virtuous cycle linking soil improvement, plant growth, and yield stability (Figures 5, 6).

Initial Phase (0–3 years): During the establishment period, soil fertility declines along a controlled trajectory similar to the fixed strategy but with smaller amplitude. The soil index drops to Soil 3 after about three years (132 ticks; Figures 5a, 6a), yet sequential fertilizer adjustments maintain robust plant growth. The fastest-growing tree species are marked to undergo a first thinning (50%), while I. paraguariensis continues to be harvested at reduced but steady increments (Figures 5b, 6b).

Transition Phase (3–5 years): The adaptive strategy produces a markedly superior soil trajectory compared with the fixed‐input case, maintaining continuous fertility recovery through optimized feedback control (Figures 6a, 7). Soil quality improves steadily, reflecting a more efficient balance between endogenous and exogenous inputs. Biomass accumulation and harvest stability (Figures 5, 8) respond with a short delay to this improvement, revealing a lagged but coherent coupling between soil restoration and aboveground productivity.

Mature Phase (5–10 years): By the end of the simulation, total harvest reaches approximately 7,000 units—showing substantially stronger annual gains than the fixed‐reposition scenario. Cumulative soil recovery approaches 25,000 units from its lowest point, following the most consistent upward trajectory across all simulations. The system exhibits remarkable uniformity in plant size, soil quality distribution, and yield variability, reflecting a stable equilibrium sustained by dynamic coupling between external fertilizers and endogenous SOM contributions (Figure 7). Locally produced SOM oscillates according to its characteristic half-life, generating subtle spatial patterns in I. paraguariensis growth; however, these variations are smoother under adaptive control, as feedback regulation dampens the amplitude of soil and growth fluctuations.

Figure 9 summarizes the comparative harvest performance across all consociated tree species, providing the quantitative foundation for model validation and demonstrating the ABM’s ability to reproduce observed empirical patterns. Across all simulations starting from Soil 1, systems without repositioning exhibited continuous soil decline, whereas systems with 20–25% repositioning transitioned to increasing soil indices. In the absence of repositioning, soil fertility steadily declines, leading to eventual system collapse. However, even minimal repositioning—approximately 20–25% of extracted nutrients during cyclical harvest and thinning routines—enables the system to evolve toward a self-sustaining, higher fertility state. In our investigation, while favorable initial soil conditions provide a head start, adaptive nutrient management remains crucial for maintaining productivity and resilience in the long term.

A quantitative comparison of the two primary management approaches is presented in the Supplementary Materials. This quantitative comparison demonstrates the adaptive strategy’s superior performance across several metrics, including higher overall productivity, increased biomass production, faster system recovery, and greater long-term soil restoration, while maintaining comparable harvest stability.

Our simulations, particularly for the ‘Soil 1’ initial condition with both fixed minimum and adaptive repositioning strategies, provide illuminating insights into the influence of individual tree species on I. paraguariensis harvest yields. As illustrated in Figure 9, distinct patterns emerged, with some tree species supporting higher I. paraguariensis harvest yields than the control (no trees), while others resulted in lower cumulative yields.

Notably, the model predicts that Toona, Cañafístola, Petiribi, and Anchico support higher harvest yields, a result consistent with observations from our real-life experimental field system (Comolli et al., 2024) over an 8-year period. This agreement strengthens the model’s validity in capturing beneficial inter-species synergies. The case of Kiri in our simulations also coincides with the field system (Comolli et al., 2025b), where it was completely cut down after three years with its biomass contributing significantly to SOM. Our current model does not simulate this specific, dramatic, management intervention for Kiri, but the harvest biomass from this fast-growing species matches this result.

For Lapacho and Grevilea our simulations predict the lowest harvest volumes, which contradict the experimental results. These results are a consequence of the slower growth rate and biomass production for these species within the model. While in the experimental system Araucaria is correlated with the highest harvest volume (Comolli et al., 2024), our simulations predict a median harvest average. These contrasts suggest possible, significant and still uncharacterized synergies. Araucaria is known to have significant synergies with I. paraguariensis (Capellari et al., 2017; Fernandez et al., 1997; Day et al., 2011; Montagnini et al., 2011, 2006), perhaps due to mycorrhiza or to enrichment of key elements, which are outside the attributes included in our simulations. Guatambu produces high volumes of leaves and likely contributes to PAR reduction more than we account for in the model.

The remaining tree species in our simulations, which also resulted in slightly lower cumulative I. paraguariensis harvest, are consistent with field experiment outcomes. This pattern is primarily attributed to the complex interplay between shade competition and localized soil improvements from tree management, which collectively balance out to these observed yield differentials.

The increasing adoption of ABMs in ecology reflects their capacity to capture feedbacks, thresholds, and emergent properties that are difficult to represent with reductionist approaches. Conventional agricultural models—whether empirical, econometric, or differential-equation based—often assume homogeneous agents and linear responses, limiting their ability to describe nonlinear, spatially explicit dynamics typical of real agroecosystems (Dupraz et al., 2019; Rahman et al., 2023). In contrast, ABMs represent individual organisms and management units as autonomous agents, allowing local interactions to generate emergent, system-level behavior (Grimm et al., 2005, 2010; Railsback and Grimm, 2012; Spies, 2017).

Within agroforestry research, relatively few ABM studies have addressed multispecies interactions over long temporal horizons. Previous work has often emphasized human–environment decision dynamics (Berger and Schreinemachers, 2011; Burgess et al., 2019) or landscape/policy applications (Torralba et al., 2018; Liu et al., 2021; Page et al., 2013; Johnson and Salemi, 2022), and many approaches simplify crop–tree biophysical coupling (Dupraz et al., 2019; Soualihou et al., 2021; Millington and Wainwright, 2016). Our model advances this literature by directly integrating species-specific growth rates, canopy shading, and soil organic matter (SOM) dynamics measured in a long-term experimental system (Comolli et al., 2024, 2025a).

Across simulations, positive feedback loops emerge when biomass generated through harvest, thinning, and pruning is returned to the soil as SOM, increasing fertility and enhancing subsequent plant growth. This biomass–SOM–fertility coupling progressively stabilizes yields and reduces dependence on external inputs. These dynamics parallel field observations of soil recovery and yield maintenance in long-term agroforestry systems (Montagnini et al., 2006; Day et al., 2011; Comolli et al., 2023; Comolli et al., 2025b; Comolli et al., 2024, 2025a) and align with theoretical descriptions of agroecosystems as complex adaptive systems governed by reinforcing feedbacks (Levin, 1998; Holling, 2001; Paterno et al., 2024).

Our agent-based model of an integrated multispecies agroforestry system demonstrates how computational approaches can elucidate complex ecological interactions in managed production landscapes. Grounding the model in empirically measured allometric parameters—including species-specific growth rates and photosynthetically active radiation (PAR)—ensures that relative plant sizes, competition intensity, and canopy interactions scale realistically. This empirical calibration enables faithful simulation of biomass production, management interventions, and their coupled effects on soil organic matter (SOM) dynamics, which underpin the model’s predictive relevance and real-world applicability.

Although soil quality is represented using dimensionless indices rather than absolute physical units, the scaling of these indices is explicitly designed to reflect observed system behavior across degraded, intermediate, and restored soil states. This abstraction allows the model to highlight emergent patterns, nonlinear responses, and threshold effects that may be difficult to isolate in purely empirical studies. As with other agent-based and process-based approaches, this level of abstraction is not a limitation but a methodological choice that facilitates the identification of key knowledge gaps, the generation of quantitative hypotheses for field testing, and the prioritization of promising system designs for experimental evaluation (Grimm et al., 2005; Railsback and Grimm, 2012).

The model captures essential spatial dynamics by representing local interactions such as shade competition, differential SOM inputs from tree management, and species-specific effects on neighboring Ilex paraguariensis plants and their harvest yields. Importantly, simulations reproduce counterintuitive field observations—such as larger I. paraguariensis individuals occurring near fast-growing tree species rather than in full-sun control rows—a pattern documented empirically (Comolli et al., 2025b) but difficult to anticipate without computational modelling.

At the same time, the simulations expose clear knowledge gaps. While the model correctly reproduces relative yield advantages for several tree species, it underestimates the positive effects of Lapacho, Grevillea, and Araucaria observed in the field. This discrepancy points to additional ecological mechanisms not yet explicitly represented—such as phosphorus mobilization, mycorrhizal associations, or facilitative microclimatic effects—and highlights concrete directions for targeted empirical and modelling follow-up.

The simulations indicate that long-term sustainability in multispecies agroforestry systems depends on the presence of critical thresholds in soil fertility and management intensity. Below these thresholds, productivity declines and self-reinforcing soil–plant feedbacks fail to establish; above them, the system converges toward a resilient equilibrium characterized by sustained biomass production and improving soil quality. In particular, the model identifies an inflection zone around a soil quality index of approximately 3,000–4,000 (corresponding to “Soil 4”–”Soil 3” conditions), below which continuous external inputs are insufficient to prevent long-term degradation.

Once this threshold is exceeded—either through prior restoration or adequate early fertilization—the system transitions to a regime dominated by endogenous SOM accumulation. In this regime, biomass inputs from harvest residues, thinning, and pruning amplify soil fertility gains, generating a positive feedback loop between productivity and soil quality. These dynamics suggest the possibility of developing empirically informed heuristics to guide site selection, initial restoration efforts, and early-stage management decisions.

Adaptive management strategies play a critical role in crossing and stabilizing these thresholds. In the simulations, feedback-regulated fertilization accelerates recovery by increasing inputs during early depletion phases and reducing them as endogenous SOM-driven fertility gains accumulate. This behavior mirrors resilience mechanisms observed in natural and managed ecosystems, where stability emerges from the balance between disturbance and recovery rates (Levin, 1998; Holling, 2001).

The model also exhibits hysteresis-like dynamics: once soil fertility falls below the viability threshold, restoring previous input levels is insufficient to recover productivity. Instead, a substantial initial investment in soil restoration is required to re-enter the self-sustaining regime. This result underscores the ecological and economic importance of preventing soil degradation beyond critical limits, as recovery costs and times increase sharply once thresholds are crossed.

Taken together, these threshold behaviors provide a conceptual bridge between empirical soil science and complex systems theory. They indicate that sustainable agroforestry systems operate within a bounded region of parameter space where soil fertility, biodiversity, and productivity co-evolve under adaptive feedback. Managing these thresholds thus emerges as a central design principle for both model-based planning and field implementation of resilient, long-lived agroecosystems.

Our simulations show that agroforestry performance is governed not only by input levels but by the timing and adaptability of management interventions. Across all evaluated metrics, adaptive repositioning—initiated at 25% of extraction and regulated by feedback—consistently outperforms fixed strategies in long-term productivity, soil restoration, and system recovery.

Adaptive management shortens the vulnerable establishment phase, enabling the system to recover from early soil depletion more than twice as fast as fixed-input strategies. Although fixed repositioning produces smoother short-term harvest trajectories, this apparent stability comes at the cost of lower cumulative productivity and slower soil recovery. In contrast, adaptive management accepts moderate temporal oscillations in exchange for faster stabilization, higher total biomass production, and substantially greater long-term yields.

These results reveal a fundamental trade-off between short-term input efficiency and long-term ecological resilience. Fixed strategies maximize harvest per unit of fertilizer, whereas adaptive strategies maximize total productivity, soil restoration, and system persistence. From a management perspective, the simulations suggest that responsiveness to system feedback—rather than static optimization—offers a more robust pathway toward sustainable agroforestry outcomes.

Importantly, the adaptive approach translates directly into actionable guidance: real-time adjustment of nutrient inputs during early depletion phases, followed by gradual reduction as endogenous SOM accumulation strengthens system fertility. Such feedback-driven strategies provide a practical framework for improving both economic viability and ecological performance in complex perennial agroecosystems.

A central emergent property of the model is the formation of self-reinforcing feedback loops linking biomass production, soil organic matter (SOM) accumulation, and plant growth. Through harvesting, thinning, and pruning, biomass is returned to the soil, increasing fertility and accelerating subsequent growth. Over time, this endogenous SOM production becomes the dominant driver of soil improvement, progressively reducing reliance on external inputs.

The inclusion of SOM decay—modelled with a finite half-life—introduces realistic temporal dynamics that prevent unrealistically smooth system trajectories. This turnover generates spatial and temporal heterogeneity in soil quality, which in turn supports resilience by providing multiple recovery pathways following disturbance. Rather than undermining stability, these controlled fluctuations contribute to long-term persistence by preventing nutrient lock-in and localized depletion, in agreement with natural systems (Pinho et al., 2012; Liu et al., 2018).

Together, these dynamics illustrate how sustainability in agroforestry systems arises not from static equilibria, but from adaptive cycling between growth, disturbance, and recovery. The model demonstrates that resilience emerges when management practices reinforce—rather than suppress—these natural feedbacks, allowing productivity and soil quality to co-evolve within a stable operational envelope.

A central implication of our results is that the focal crop, Ilex paraguariensis, does not merely respond to its associated tree species but actively reshapes their collective productive environment through feedbacks mediated by harvest, pruning, and soil organic matter (SOM) dynamics. Tree species modify light availability, microclimate, and nutrient inputs to the crop, while the crop’s biomass extraction and residue return—together with tree biomass converted to SOM through crop-driven management—reciprocally reorganize soil fertility gradients that govern tree regrowth and competitive balance. Productivity in this system therefore emerges as a self-reinforcing dynamic, a network property, rather than as an additive contribution of independent species. This explicitly demonstrates that agroforestry yield cannot be attributed to single-species performance alone, but instead to the structure of inter-species coupling and management-mediated feedback. Agroforestry systems function as mutualistic collectives rather than hierarchical arrangements. The model provides a computational framework for exploring how these reciprocal dependencies scale up to generate system-level properties like resilience, stability, and sustained productivity. Future empirical work should explicitly test these bidirectional couplings through controlled variation of tree and I. paraguariensis management regimes. The model situates agroforestry within the broader field of complex systems science, bridging empirical ecology, systems modelling, and sustainability policy.

While the present ABM captures essential biophysical feedbacks and emergent properties of multispecies agroforestry systems, it necessarily simplifies several ecological processes. The model does not yet explicitly represent microbial communities, mycorrhizal networks, faunal interactions, or detailed biogeochemical fluxes, all of which are known to influence nutrient cycling and ecosystem resilience. Future extensions could incorporate climate variability, pest–predator dynamics, and socio-economic feedbacks to enhance realism and broaden applicability across agroecological contexts. Importantly, the current formulation is intentionally minimal and empirically grounded, allowing it to expose causal structure and knowledge gaps that can guide targeted experimental and modelling efforts. Beyond its immediate relevance to yerba mate–based agroforestry, the framework developed here can be adapted to other multispecies agricultural systems and tree planting projects (Francini et al., 2024; Schwaab et al., 2021; Paterno et al., 2024; Johnson and Salemi, 2022; Spies, 2017).