Hi-sAFe: A 3D Agroforestry Model for Integrating Dynamic Tree–Crop Interactions

: Agroforestry, the intentional integration of trees with crops and / or livestock, can lead to multiple economic and ecological beneﬁts compared to trees and crops / livestock grown separately. Field experimentation has been the primary approach to understanding the tree–crop interactions inherent in agroforestry. However, the number of ﬁeld experiments has been limited by slow tree maturation and di ﬃ culty in obtaining consistent funding. Models have the potential to overcome these hurdles and rapidly advance understanding of agroforestry systems. Hi-sAFe is a mechanistic, biophysical model designed to explore the interactions within agroforestry systems that mix trees with crops. The model couples the pre-existing STICS crop model to a new tree model that includes several plasticity mechanisms responsive to tree–tree and tree–crop competition for light, water, and nitrogen. Monoculture crop and tree systems can also be simulated, enabling calculation of the land equivalent ratio. The model’s 3D and spatially explicit form is key for accurately representing many competition and facilitation processes. Hi-sAFe is a novel tool for exploring agroforestry designs (e.g., tree spacing, crop type, tree row orientation), management strategies (e.g., thinning, branch pruning, root pruning, fertilization, irrigation), and responses to environmental variation (e.g., latitude, climate change, soil depth, soil structure and fertility, ﬂuctuating water table). By improving our understanding of the complex interactions within agroforestry systems, Hi-sAFe can ultimately facilitate adoption of agroforestry as a sustainable land-use practice.


Introduction
Agroforestry, the intentional integration of trees with crops or livestock, has gained momentum in temperate regions over the last few decades [1][2][3]. Temperate agroforestry research, however, has lagged behind work in the tropics and subtropics [4,5]. At the core of agroforestry research are the interactions between trees and crops, which drive differences from monoculture cropping systems and pure forestry stands. These interactions can lead to both positive and negative effects when mixing trees and crops in agroforestry. In many environments, interaction between trees and crops can result in overyielding in agroforestry compared to tree and crop monocultures [6][7][8]. Agroforestry system design can exploit these benefits to create economic advantages [9][10][11], but requires specific knowledge on the complex tree-crop interactions.
Field experimentation has been the primary approach to understanding the complex dynamics within agroforestry systems. Well-designed, long-term experiments can produce data for many years, but slow tree maturation, difficulty in obtaining consistent funding, and other logistical constraints have hindered success [12]. The scope of temperate agroforestry research has been constrained by the relatively few long-term field experiments. For example, 30% of temperate publications on alley cropping (i.e., silvoarable) field experiments over the last three decades originated from just three sites [5]: the Guelph Agroforestry Research Station (Guelph, Canada; est. 1988) [13], Domaine de Restinclières (Montpellier, France; est. 1995) [14], and the University of Missouri Greenley Memorial Research Center (Novelty, MO, USA; est. 1997) [15].
Models of agroforestry systems have the potential to advance ecological understanding while providing greater guidance for future experimentation [16]. However, the complexity of interactions in agroforestry systems in both space and time makes the development of models difficult. Key goals for agroforestry models include the capacity to simulate: • systems across pedoclimatic environments and management regimes; • above-and belowground interactions for light, water, and nutrients; • the multiple possible yields, including food, fiber, and fuel; and • ecosystem services such as the capture of excess nutrients, soil erosion, and carbon sequestration [17].
Process-based rather than empirical models are the ideal choice for agroforestry systems as they better mimic the complex dynamics of tree-crop interactions and allow predictions outside the range of data available for model parameterization [18][19][20]. The availability of process-based agroforestry models is limited, however, despite a rich array of such models for both forestry and agronomy [16,17].
Simple agroforestry models were originally modified from existing crop models. For example, CROPGRO [21] and STICS [22] have been used to model agroforestry systems by simply reducing the amount of light available to crops [23,24]. Similarly, CROPGRO and EPIC [25] were used to evaluate the effect of shelterbelts on crops by modifying crop exposure to wind and radiation [26,27]. The WIMISA model also integrates belowground competition for water between trees and crops [28]. None of these simple models, however, predict tree growth nor model crop productivity for multiple seasons.
Some process-based models have specifically been developed to handle agroforestry systems. Yield-SAFE, for example, is a parameter-sparse, one-dimensional (1D) biophysical model developed to simulate the productivity of an agroforestry system for an entire tree rotation [29]. The model considers light and water interactions and is coupled to a bio-economic model for evaluating profitability [6,30]. Huth et al. [31] proposed an adaptation of the APSIM crop model to simulate growth of a hedge and its competition for water and light with an adjacent crop. This model is two-dimensional (2D), accounting for interactions as a function of distance from the hedge, but no interactions within the hedge. The HyPAR model [32], a coupling of the Hybrid forest model [33] and PARCH crop model [34], has been used to predict the productivity of agroforestry systems along a wide gradient of aridity [35]. The most commonly and successfully used agroforestry model is WaNuLCAS [36], which is the only model capable of integrating light, water, and nitrogen competition for an entire system rotation. WaNuLCAS has been used successfully in a range of tropical agroforestry systems [37][38][39][40], but was not designed for temperate systems. Application of all agroforestry models to date has remained limited due to inadequate flexibility, restricted complexity in simulating interactions, or difficult parameterization requirements [17].
Here, we introduce Hi-sAFe, a three-dimensional (3D), process-based, biophysical model that integrates tree-crop interactions in agroforestry systems. Hi-sAFe has been under development since 2002 via the Silvoarable Agroforestry for Europe (SAFE) project [41] and was partially described by Talbot [42]. The model attempts to overcome the gaps and weaknesses of existing agroforestry models by capturing spatial and temporal heterogeneity. Two main hypotheses directed Hi-sAFe development: • Productivity of an agroforestry system depends on the acquisition of heterogeneously distributed resources by trees and crops [43]. • Tree-crop interactions are, in part, governed by aboveground, belowground, and phenological plasticity [44].
To explore these hypotheses, the specific objectives of Hi-sAFe development were to create a model that could simulate: • three-dimensional tree-crop interactions for light, water, and nitrogen; • plastic aboveground and belowground tree architecture responsiveness to resource availability; • the full lifetime of the system, from tree planting to harvest, on a daily time-step; and • the principal agroforestry design and management strategies, such as branch pruning, tree thinning, root pruning, and the incorporation of an uncropped area around each tree or strip along the tree row.
This paper synthesizes the approach of Hi-sAFe and provides a complete model description. We also discuss model applications, limitations, and implementation.

Model Specification
Hi-sAFe simulates a 3D agroforestry environment by coupling the existing STICS (version 8) crop model [22,45] with a new tree model (sAFe-Tree) on a daily time-step ( Figure 1). The integration of STICS and sAFe-Tree occurs via three tree-crop interaction modules that govern the dynamics of light, water, and nitrogen via simple, well-established equations. Hi-sAFe is written in Java and runs on the Capsis platform [46]. Here, we introduce Hi-sAFe, a three-dimensional (3D), process-based, biophysical model that integrates tree-crop interactions in agroforestry systems. Hi-sAFe has been under development since 2002 via the Silvoarable Agroforestry for Europe (SAFE) project [41] and was partially described by Talbot [42]. The model attempts to overcome the gaps and weaknesses of existing agroforestry models by capturing spatial and temporal heterogeneity. Two main hypotheses directed Hi-sAFe development: • Productivity of an agroforestry system depends on the acquisition of heterogeneously distributed resources by trees and crops [43].
To explore these hypotheses, the specific objectives of Hi-sAFe development were to create a model that could simulate: • three-dimensional tree-crop interactions for light, water, and nitrogen; • plastic aboveground and belowground tree architecture responsiveness to resource availability; • the full lifetime of the system, from tree planting to harvest, on a daily time-step; and • the principal agroforestry design and management strategies, such as branch pruning, tree thinning, root pruning, and the incorporation of an uncropped area around each tree or strip along the tree row.
This paper synthesizes the approach of Hi-sAFe and provides a complete model description. We also discuss model applications, limitations, and implementation.

Model Specification
Hi-sAFe simulates a 3D agroforestry environment by coupling the existing STICS (version 8) crop model [22,45] with a new tree model (sAFe-Tree) on a daily time-step ( Figure 1). The integration of STICS and sAFe-Tree occurs via three tree-crop interaction modules that govern the dynamics of light, water, and nitrogen via simple, well-established equations. Hi-sAFe is written in Java and runs on the Capsis platform [46].

Hi-sAFe Scene
The model simulates a rectangular scene (Figure 2), which can be replicated in all directions via periodic boundary conditions (PBC; see Figure 3). The scene is divided into a grid of square "cells", each of which grows a homogenous crop and up to one tree of any number of species. The positions of trees define two areas, the cropped area (usually an "alley" if the scene is replicated), and the uncropped area close to the trees, each of which may grow one of several plant species or be bare soil. Cells are vertically divided into discrete "voxels" to capture belowground soil structure, with soil properties (root density and water and nitrogen content) homogenous within a voxel. Up to five "layers" of physical soil properties can be applied, each homogenous within the layer. Soil properties for each layer include sand, silt, clay, limestone, organic matter, and stone content. Cell, voxel, and layer dimensions can be modified to match field conditions and manage the spatial resolution of model predictions. The sky is represented by a hemispherical vault discretized into sectors defined by elevation and azimuth intervals. Hi-sAFe scene description parameters are shown in Table 1.

Hi-sAFe Scene
The model simulates a rectangular scene (Figure 2), which can be replicated in all directions via periodic boundary conditions (PBC; see Figure 3). The scene is divided into a grid of square "cells", each of which grows a homogenous crop and up to one tree of any number of species. The positions of trees define two areas, the cropped area (usually an "alley" if the scene is replicated), and the uncropped area close to the trees, each of which may grow one of several plant species or be bare soil. Cells are vertically divided into discrete "voxels" to capture belowground soil structure, with soil properties (root density and water and nitrogen content) homogenous within a voxel. Up to five "layers" of physical soil properties can be applied, each homogenous within the layer. Soil properties for each layer include sand, silt, clay, limestone, organic matter, and stone content. Cell, voxel, and layer dimensions can be modified to match field conditions and manage the spatial resolution of model predictions. The sky is represented by a hemispherical vault discretized into sectors defined by elevation and azimuth intervals. Hi-sAFe scene description parameters are shown in Table 1. The simulated Hi-sAFe scene reduces a landscape into its simplest replicable unit and is discretized into a 3D, spatially explicit voxel grid. Periodic boundary conditions (PBC) allow for the scene to be repeated infinitely in all directions. When no model constraints limit the range of a parameter, "-" is used. Square brackets indicate that the interval endpoint is included in the allowed values, whereas parentheses indicate that the interval endpoint is excluded from the allowed values. The simulated Hi-sAFe scene reduces a landscape into its simplest replicable unit and is discretized into a 3D, spatially explicit voxel grid. Periodic boundary conditions (PBC) allow for the scene to be repeated infinitely in all directions.

STICS Crop Model in Hi-sAFe
The STICS crop model simulates crop growth, crop management interventions (e.g., fertilization, tillage, irrigation), 1D vertical soil fluxes of water and nitrogen, and soil organic matter processes. Horizontal soil water and nutrient fluxes are not modeled. To integrate STICS within a 3D heterogeneous environment, Hi-sAFe runs an instance of STICS for each cell in the simulated scene. Aggregation and disaggregation algorithms are used for communication between the 1-cm thick mini-layers of STICS and the~20-cm thick voxels of Hi-sAFe. The required daily climate inputs of Hi-sAFe are identical to those of STICS: minimum and maximum air temperature, minimum and maximum relative air humidity, global radiation, precipitation, and wind speed. In addition, a daily water table depth can be provided to simulate a fluctuating water table. Generic parameters passed to STICS are not discussed here.

Water and Nitrogen Cycles
The water cycle in Hi-sAFe includes inputs to the scene via rain, snow, irrigation, and a dynamic water table. The water table is simulated by maintaining voxels below a given depth saturated with water and can provide water to the scene by (1) adding water to a voxel when bringing it to saturation and (2) directly supplying roots in a saturated voxel with an inexhaustible supply of water. Water exports from the scene include soil surface runoff, soil evaporation, natural and artificial drainage, and crop and tree transpiration.
The nitrogen cycle includes inputs to the scene via fertilization, irrigation, deposition, fixation by crops, and the water table. Nitrogen fixation by trees is not modeled. The water table maintains a fixed nitrate concentration and can provide nitrogen to the scene by (1) adding nitrate to a voxel that has a nitrate content lower than the water table and (2) directly supplying roots in a saturated voxel with an inexhaustible supply of nitrate. Nitrogen exports from the scene include crop and tree harvests, leaching via natural and artificial drainage, nitrification, denitrification, and volatilization of fertilizer. Management of humus, crop residues (mineralization, humification, etc.), and the fate of nitrogen fertilizers are simulated by STICS.

sAFe-Tree Model
The sAFe-Tree model, inspired by the HyPAR model [32], simulates key aspects of tree growth in response to light, water, and nitrogen resources: phenology, carbon assimilation, carbon allocation among tree compartments, aboveground growth and allometry, belowground growth, and optional management interventions such as branch and root pruning. An overview of the sAFe-Tree model is described below, with all model parameters defined in Table 2. Threshold for frost mortality of leaves Duration of leaf fall days [1,365] Light interception µ Correction parameter to account for leaf clumping. A clumping coefficient below 1, equal to 1, or above 1 indicates a clumped, random, or regular distribution of leaves inside the crown, respectively. - Wood area density for light interception by tree branches Carbon allocation a NSC Threshold imbalance above which remobilization of C NSC is triggered Limits maximum daily ∆C NSC as a fraction of C leaves -(0, 1] α NSC * Target C NSC as a fraction of the tree woody C pool -(0, 1] LFR* 0 Initial LFR* at tree planting kg kg −1 (0, 1) Sensitivity of LFR* to N satisfaction when there is no stress Maximum allowed LFR* kg kg −1 (0, 1)

Growth and senescence SLM
Leaf dry mass per unit leaf area kg m −2 (0, -) θ leaves Leaf C content g C g −1 dry biomass (0, 1) θ wood C content of all compartments except leaves g C g −1 dry biomass (0, 1) ρ wood Wood density of branches, stems, stump, and coarse roots kg m Maximum fraction of intercepted rain routed to stem flow -

Tree Architecture
Aboveground tree architecture is composed of a trunk and stump (truncated cone) surmounted by a homogeneous canopy containing leaves and branches ( Figure 4). Belowground tree compartments include coarse roots, represented as 3D structural elements, and fine roots, quantified as a length present in each voxel ( Figure 5). Trees are initialized with a specified H, H p , r x , r y , and leaf-to-fine-root ratio.      Tree coarse root topology simulated by Hi-sAFe after 12 years of growth in an alley cropping system. Tree spacing between rows is 13 m and spacing within rows is 7 m. The simulated tree line extends into the page. The vertical axis is stretched by three times for visualization clarity. Illustrated coarse root diameter reflects actual coarse root diameter. Coarse root color is proportional to fine root density in that voxel, with darker colors indicating more fine roots. Roots primarily remain below the crop rooting depth, with surface roots only occurring within the untilled tree line. Some roots can also be seen growing back into the tilled alley crop area from the untilled tree line and from roots below the tillage depth.

Phenology
Hi-sAFe was designed for deciduous trees typical of temperate climates. Trees are dormant until day D 0 , after which bud burst is triggered following the accumulation of S T degree-days. Leaf expansion continues for D LE days, with daily dynamics managed by the carbon allocation algorithms presented below. Leaf fall begins at the onset of physiological senescence D V days after bud burst or if the average daily temperature drops below T LF . Leaf fall occurs over D LF days following a decreasing sigmoid function. The tree reenters dormancy once leaf fall is complete.

Light Interception
Light interception by trees is calculated via a ray tracing algorithm adapted from Courbaud et al. [48]. Ray tracing models are common in other 3D biophysical models, such as MIR [49], tRAYci [50], and SORTIE [51]. The position of the sun is calculated every half hour as a function of latitude, date, cardinal orientation of the scene, and slope of the scene. Direct radiation from each period is attributed to the closest sky vault sector. Diffuse radiation is calculated following Moon and Spencer [52] and divided across all sectors proportionally to their area. Daily radiation in each sector is summed across all periods in the day, and the trajectory of light rays between each sky sector and scene cell is calculated. The fraction of incident light in a ray transmitted through a tree crown (L trans /L 0 ) is calculated as: where L is the length of the ray path through the crown, µ is a leaf clumping coefficient [42], σ is the leaf absorbance for the considered wavelengths, LAD is the leaf area density within the crown, and WAD is the wood area density within the crown [53]. The model assumes that transmittance and reflection of rays by trunks and branches are zero, but branches do not shade leaves of the same tree. Therefore, the fraction of incident light absorbed by tree leaves (L leaves /L 0 ) is calculated as: Ray interception is integrated over the scene to obtain the radiation and photosynthetically active radiation (PAR) intercepted by individual trees and transmitted to each cell in the scene. With PBC enabled, shading is transmitted to the opposite side of the scene (Figure 3b).

Carbon Assimilation
Net carbon assimilation by the tree on day j (A j ) is modeled as an empirical function of light-use efficiency (LUE j ) and the intercepted PAR (PAR i,j ) [54]: LUE j is calculated as a function of a maximum potential LUE for the species (LUE max ) reduced by a water stress index (W stress,j ), a nitrogen stress index (N stress,j ), and a leaf age effect (LA j ): where ω LUE is the sensitivity of LUE to water stress and η LUE is the sensitivity of LUE to nitrogen stress. Water stress is calculated as the ratio of water uptake (W uptake,j ) to water demand (W demand,j ), and nitrogen stress is calculated as the ratio of total tree nitrogen content (N total,j ) to optimum tree nitrogen content (N optimum,j ). Both W stress,j and N stress,j , therefore, range from zero (complete stress) to one (no stress). The sensitivity factors ω LUE and η LUE can take on any non-negative value, with zero representing no responsiveness to the respective stress index. The leaf age effect is calculated as: where λ is a leaf senescence time constant, D BB is the date of bud burst, and τ is the leaf age at which maximum LUE occurs [55].

Carbon Allocation
The tree contains four aboveground carbon pools (leaves-C leaves , branches-C branches , stem-C stem , and stump-C stump ; collectively C AG ), two belowground carbon pools (coarse roots-C cr and fine roots-C fr ; collectively C BG ), and a labile pool of nonstructural carbon (C NSC ). Carbon allocated to tree growth on day j (C G,j ) is equal to A j plus ∆C NSC,j . Labile carbon is remobilized to increase C G,j if C leaves,j < (1 − a NSC ) * C leaves * (6) where C leaves * is the target leaf carbon pool size calculated via Equation (A5) ( Table 3) and a NSC is an imbalance threshold. To prevent the entire labile pool from being remobilized on one day, ∆C NSC,j is calculated as: where b NSC and c NSC are limiting fractions of existing carbon pools. In contrast, replenishment of C NSC by C G,j occurs if C NSC,j < C braches,j + C stem,j + C stump,j + C cr,j * α NSC * where α NSC * is a target proportion of woody carbon pools. Allocation of C G,j between aboveground and belowground compartments is governed by a target ratio between leaf and fine root carbon pools (LFR*), based on the concept that the compartment limiting primary productivity must be allocated more carbon to better exploit the limiting resource (functional equilibrium): Each tree is initialized with LFR 0 *, and then LFR* changes daily by a maximum amount of ∆LFR * max . In the absence of water and nitrogen stress, the tree favors aboveground growth, and LFR* increases following a coefficient δ LFR , the nitrogen saturation of the tree (N sat,j ; Equation (17)), and a sensitivity factor (ε LFR ) to N sat : Conversely, if water or nitrogen stress occur, a reduction in LFR* is calculated as a function of the stresses and stress sensitivities, in a manner analogous to that for LUE j (Equation (4)): where ω LFR is the sensitivity of LFR* to water stress and η LFR is the sensitivity of LFR* to nitrogen stress. In all cases, LFR* is constrained by an upper (LFR * max ) and lower (LFR * min ) limit. The fraction of C G,j allocated to aboveground compartments (F AG,j ) is then calculated via: following the functional equilibrium as described by Poorter and Nagel [56].

Growth and Allometry
All carbon allocated to aboveground growth (C G,j × F AG,j ) is initially allocated to C stem . Using this temporary pool, new target aboveground tree architecture is calculated via seven allometric relationships (Equations (A1)-(A7); Table 3) and fixed values for specific leaf mass (SLM), leaf carbon content (θ leaves ), wood carbon content (θ wood ), and wood density (ρ wood ). The carbon sink strength for each compartment is then calculated as the difference in the pool size between the target tree and current tree. Actual carbon allocation among compartments is proportional to their individual carbon sink strength. Finally, aboveground architecture dimensions are updated using the actual carbon allocation to each compartment, with horizontal crown radii stopping once the crown touches another tree crown. Table 3. Allometric relationships governing tree aboveground architecture in sAFe-Tree. Allometric parameters are not included in Table 2.
V stem : stem volume; V crown : crown volume; V branches : branch volume; V stump : stump volume; A crown : area covered by the projection of the canopy onto the soil surface; A leaves : leaf area.
Carbon allocation among belowground compartments is governed by an algorithm based on both local (voxel-level) and global (tree-level) controls of root growth, with fine root development driving coarse root architecture [47]. Carbon allocation to fine root proliferation within a rooted voxel and colonization of adjacent voxels depends on (1) the success of water and nitrogen uptake in the voxel during the previous day, relative to overall supply and demand at the tree level, and (2) the topological distance from the voxel to the stem base. The seven parameters of the fine root growth algorithm are  Table 2, and algorithm details can be found in Mulia et al. [47]. Coarse root growth is then determined by "pipe-stem" model principles and a fixed ratio of coarse root cross-sectional area to fine root length (γ). With PBC enabled, roots can grow out of and re-enter the opposite side of the scene (Figure 3c). Fine roots have a fixed specific root length (SRL), and both root compartments have a fixed carbon content (θ wood ).
Fine root senescence follows an exponential decay with mean lifespan τ fr . Anoxia can occur within voxels either when the water table is above voxel depth or when slow infiltration temporarily leads to a perched water table. In anoxic conditions, fine root mean lifespan decreases to τ fr,anoxia . Anoxia can also cause coarse root death, which occurs after τ cr days of anoxic conditions. When a segment of coarse root dies, all downstream coarse and fine roots also die. All root litter enters the soil with their carbon and nitrogen concentrations at senescence.

Water Demand and Rain Interception
The daily evaporative demand of a tree is calculated as: where Γ is a scaling coefficient, f tree,j is the proportion of incoming radiation intercepted by the tree, and PET 0 is the potential evapotranspiration calculated from meteorological inputs. Tree water uptake (WU tree,j ) is calculated via the tree-crop interaction algorithm described below. Rain stored on tree leaves (P stored,j ), which accumulates from incident rain (P incident,j ) that is intercepted by tree crowns, can reduce ED tree,j when it evaporates. Intercepted rain (P intercepted,j ) is calculated as: where ω is leaf wettability and LAI j is the leaf area index of the tree. Intercepted rain is either transferred to stem flow or evaporated off leaves, thereby also reducing the rain that reaches the cells directly below the tree crown. The amount of intercepted rain directed to stem flow (SF j ) is calculated as: where SF max is the maximum fraction of rain that leaves as stem flow and c SF is a scaling coefficient.

Nitrogen Demand and Allocation
Tree nitrogen demand (N demand,j ) is calculated to attempt to maintain (1) optimum nitrogen/carbon ratios (η c ) for each structural compartment (c) following carbon allocation and (2) a nitrogen reserve: where N total,j is the total nitrogen content in all tree compartments and α N is a scaling coefficient defining the target nitrogen reserve. Nitrogen saturation, an index used to influence changes in LFR*, is calculated as: where N optim is the total of the optimum nitrogen contents in structural compartments and β N is a coefficient defining "luxury" nitrogen content. Tree nitrogen uptake (NU tree,j ) is calculated via the tree-crop interaction algorithm described below. Nitrogen allocation among structural compartments is proportional to their individual nitrogen sink strength. The nitrogen sink strength for each compartment (N sink,c,j ) is calculated as the difference between the optimum and actual (N c,j ) pool size: Total nitrogen allocated to structural compartments is then calculated as: Any excess N uptake,j not allocated to structural compartments is allocated to N labile . To conserve nitrogen, fixed fractions of leaf (κ leaves ) and fine root (κ fr ) nitrogen are remobilized to N labile prior to senescence.

Tree Management Interventions
Basic agronomic management interventions, such as fertilization, irrigation, tillage, and mulching, are managed by STICS for each cell in the Hi-sAFe scene. Three additional tree management interventions are included in Hi-sAFe. Tree thinning kills a tree and removes all aboveground biomass. Branch pruning can (1) modify r x , r y , and r z based on pruning to an H p specified as a fraction of H or a maximum allowed H p (r x and r y are reduced proportionally to the reduction in r z ) and/or (2) maintain H p and crown radii while reducing branch and leaf biomass by a specified proportion ( Figure 6). The first branch pruning approach is primarily used on timber trees to maintain long, branch-free stems, whereas the latter is primarily used in fruiting trees to enhance airflow and fruit quality. Finally, root pruning is specified by a pruning depth and distance from the tree line. For any coarse roots that are cut at this location, all downstream coarse and fine roots are killed ( Figure 6).  Fine root density is proportional to voxel shading, with darker colors indicating more fine roots. Branch pruning to a height Hp reduces vertical and horizontal crown radii by the same proportion. A reduction in WAD (and consequently LAD) can also be specified. Root pruning occurs along equidistant, parallel lines that straddle the tree (zig-zag lines; into the page). Coarse roots that are cut by root pruning (dashed lines) are killed, along with all downstream coarse and fine roots (hatched voxels). It is possible for vertically growing coarse roots to avoid root pruning and maintain roots above the pruning depth, as shown on the left side of the illustrated scene. LAD: leaf area density within the crown; WAD: wood area density within the crown.

Water Competition and Uptake
The algorithm for tree-crop competition for water uptake extends the approach of de Willigen et al. [57] to include any number of plants with overlapping root systems. This algorithm requires seven parameters for each plant and provides a global solution to water uptake without assuming an advantage by one species, taking full advantage of the dynamic tree root topology in Hi-sAFe. Full details on the algorithm are available on the Hi-sAFe website [58], but a summary is provided below. The algorithm is similar to one in WaNuLCas [36], but is extended to be 3D.
Competition for water involves the response of each plant to the available water supply in their Fine root density is proportional to voxel shading, with darker colors indicating more fine roots. Branch pruning to a height Hp reduces vertical and horizontal crown radii by the same proportion. A reduction in WAD (and consequently LAD) can also be specified. Root pruning occurs along equidistant, parallel lines that straddle the tree (zig-zag lines; into the page). Coarse roots that are cut by root pruning (dashed lines) are killed, along with all downstream coarse and fine roots (hatched voxels). It is possible for vertically growing coarse roots to avoid root pruning and maintain roots above the pruning depth, as shown on the left side of the illustrated scene. LAD: leaf area density within the crown; WAD: wood area density within the crown.

Water Competition and Uptake
The algorithm for tree-crop competition for water uptake extends the approach of de Willigen et al. [57] to include any number of plants with overlapping root systems. This algorithm requires seven parameters for each plant and provides a global solution to water uptake without assuming an advantage by one species, taking full advantage of the dynamic tree root topology in Hi-sAFe. Full details on the algorithm are available on the Hi-sAFe website [58], but a summary is provided below. The algorithm is similar to one in WaNuLCas [36], but is extended to be 3D.
Competition for water involves the response of each plant to the available water supply in their respective root zones. Plant water uptake is driven by the transpiration demand and the potential to meet demand using "minimum energy" from the voxels where the plant has roots. Each plant adjusts its water demand to the available supply by lowering their water potential and closing stomata.
The steady-state solution for nutrient transport towards regularly distributed roots in a unit soil volume also applies to water if the matrix flux potential, rather than soil water content or soil water potential, is used to reflect water availability [59]. The matrix flux potential, Φ, is the integral of the unsaturated hydraulic conductivity and can be used to predict maximum flow rates through soil [60]. Conductivity depends on soil water content, transport demand, volume of soil per unit root length, and root diameter.
The algorithm begins by estimating the plant water potential required to meet the transpiration demand, given the "perceived" water potential of the rooted voxels. In the weighting of different soil layers, a drought signaling effect can be represented, whereby roots in dry layers can signal partial stomatal closure [61,62]. The initial estimate of plant water potential at the stem base starts from this soil potential and is decreased by root entry and root longitudinal transport resistances. Actual demand is reduced based on a sigmoidal function (see [57]) due to stomatal adjustments at this plant's water potential. The reduced water demand then adjusts transport resistances, and the rhizosphere potentials are differentiated over the voxels based on differences in longitudinal resistance. In some voxels, a plant may have a rhizosphere potential equal to or higher than the soil potential and, consequently, will not be able to extract water. Roots from different plants can have different rhizosphere potentials in the same voxel.
The total difference in matrix flux potential Φ Soil − Φ i between the bulk soil and the rhizosphere of the ith most negative plant of n plants rooted in the voxel can be partitioned into a series of ranges, A gradually decreasing number of plants operate in these ranges until the last range, which is only accessible to the plant with the most negative rhizosphere potential. The rank order of plants can differ among voxels due to differences in plant location and longitudinal resistances. The potential uptake U of each plant i from each voxel j is then summed over each range k: where D i,j is the root length density of plant i in voxel j, V j is the volume of voxel j, and ρ i,j is calculated as:

Nitrogen Competition and Uptake
The algorithm for tree-crop competition for nitrogen uptake in Hi-sAFe was adapted from the algorithm for zero-sink nutrient uptake by a plant root used in the WaNuLCAS model [36]. In the WaNuLCAS algorithm, plant demand is calculated as the nitrogen deficit in all plant organs. Nitrogen transport in soil includes processes related to root structural and physiological characteristics, as well as the impact of the soil water content on diffusion. Nitrogen uptake for each plant in a voxel is determined by summing the root length density and potential zero-sink nutrient supply of all plants present [60] and then sharing the potential uptake among species proportional to their root length density. While absorption mechanisms differ for ammonium and nitrate, plant uptake is modeled using their summed potential [59]. In adapting the WaNuLCAS algorithm, three more considerations were added:

1.
A weighted mean root diameter was used to account for differences in root diameter among plant species; 2.
To account for differences in current demand that may make the zero-sink assumption an overestimation for some plants, current demand per unit root length at the plant level was used as an additional weighing factor; 3.
To prevent plant uptake in a mixture from being more than it would be in a monoculture, a series of constraints were added to ensure that including roots of a nondemanding plant does not increase uptake by others.
The algorithm ensures two levels of compensation. At the plant level, if uptake opportunity in one voxel increases (via increased nitrogen concentration, increased soil water content facilitating diffusion, increased rooting), uptake from other voxels will only decrease if total plant demand has been met. At the community level, if one plant increases root length density in a single voxel, uptake by other plants may adjust in all other voxels as well.
Overall, the algorithm allows a flexible response of each plant to adjust uptake in each voxel based on nitrogen availability across rooted voxels. This can result in complex, nonmonotonic responses by each plant to changes in supply anywhere in the soil profile. In addition to competition among plants, a positive impact of resource sharing on increased overall nutrient capture can occur. Full details on the nitrogen competition Hi-sAFe algorithm are available on the Hi-sAFe website [63] and are the 3D extension of the WaNulCAS [36] algorithm.

Model Approach
The approach used in developing Hi-sAFe was to implement the relevant competition processes at their appropriate scales. Consequently, different levels of complexity and disaggregation were chosen for different components and processes. Calculation of light competition at the hourly scale and water and nitrogen competition at the daily scale proved effective in simulating the emergent properties of agroforestry systems at multiple time scales (e.g., season, year, decades). The Hi-sAFe scene was also designed to accommodate any level of system complexity, from the most common alley cropping systems containing a single tree species and single crop species [5] to innovative multispecies and multistrata systems [12,64]. Independent representation of crops on each cell also allows for precision crop management.
For some components, simple representation was deemed sufficient. For example, tree canopies are approximated as simple ellipsoids rather than with explicit representation of individual branches and leaves. Similarly, tree carbon capture uses a light-use efficiency approach rather than a mechanistic photosynthetic model. In contrast, more complex algorithms were deemed necessary for other model components. To adequately capture the spatial complexity of belowground interactions and simulate the impact of root pruning, structurally explicit tree root systems are represented at the voxel scale. Furthermore, detailed mechanisms of a fluctuating water table on water and nitrogen availability and root mortality were included because field studies have illustrated a strong impact of a water table on tree-crop interactions [65].

Model Applications
Given the paucity of long-term, temperate agroforestry field experiments, perhaps the most powerful application of Hi-sAFe is the extrapolation of results from existing experiments to longer time scales, new agroforestry designs, management strategies, and pedoclimatic environments. Existing data can be used to parameterize and calibrate the model, and then different soil specifications or climate series can easily be applied. Comprehensive model outputs on carbon, light, water, and nitrogen pools and fluxes can facilitate a wide range of analyses on complex tree-crop interactions at daily, annual, and rotation time scales (Figure 7). Exploring future scenarios under climate change is a particularly intriguing application. Agroforestry has great potential as a tool for climate change mitigation and adaptation [64,66,67]. Nevertheless, the complexities of agroforestry systems have restricted field experiments to primarily exploring individual adaptation mechanisms on crop growth (e.g., [68][69][70]). Hi-sAFe can integrate a range of potential climate change impacts, adaptation mechanisms, and their interactions with both tree and crop growth.
Overyielding in agroforestry systems relative to tree and crop monocultures is central to their Exploring future scenarios under climate change is a particularly intriguing application. Agroforestry has great potential as a tool for climate change mitigation and adaptation [64,66,67]. Nevertheless, the complexities of agroforestry systems have restricted field experiments to primarily exploring individual adaptation mechanisms on crop growth (e.g., [68][69][70]). Hi-sAFe can integrate a range of potential climate change impacts, adaptation mechanisms, and their interactions with both tree and crop growth.
Overyielding in agroforestry systems relative to tree and crop monocultures is central to their interest and potential benefits [6,8,16,44]. Potential mechanisms of overyielding include niche complementarity (i.e., interspecific differences in resource utilization), facilitative interactions among species (e.g., species utilize nitrogen fixed by legumes), and reductions in negative plant-soil feedbacks [71][72][73]. While the specific mechanisms for overyielding have been explored in crop-only [74,75] or tree-only [76][77][78][79] systems, mechanisms that drive overyielding in agroforestry have been difficult to disentangle as very few experimental sites have maintained tree and crop monocultures adjacent to agroforestry plots for the complete tree rotation [5,80]. Moreover, mechanisms of overyielding can be highly dynamic both within and among years. For example, artificial shading can have no effect on maize productivity in a dry year, but a strong negative effect in a wet year [81]. Similarly, Gea-Izquierdo et al. [82] showed that the effect of holm oaks on the growth of pastures in Mediterranean silvopastoral systems can change from positive to neutral and negative within a growing season. These dynamics are especially true in agroforestry systems, which are in perpetual transition between a state where trees have negligible influence on crops and a state where trees dominate the plot.
Process-based models can facilitate studies on the mechanisms of overyielding by joining experimental and conceptual information on tree-crop interactions. For example, Hi-sAFe can provide insights on key interactions and limiting processes via intermediate variables that are difficult to measure in the field. Comprehensive model outputs can facilitate a wide range of analyses on overyielding, as well as calculation of the land equivalent ratio (LER; [83]) using other metrics such as light capture, water uptake, and nitrogen uptake ( Figure 8). Furthermore, since Hi-sAFe can simulate tree and crop monocultures in addition to agroforestry, it provides a powerful opportunity to leverage existing field experiments that do not have all three systems in place. An existing agroforestry plot can be used for parameterization and calibration, and then the model can provide yield estimates for the respective monoculture systems.
Even when robust calibration and validation data are lacking, process-based models can be used to (1) identify and prioritize knowledge gaps, (2) generate quantitative hypotheses for testing in the field, and (3) select the most promising system designs to include in field experiments. To these ends, Hi-sAFe can be used to explore agroforestry designs (e.g., tree spacing, crop type, tree row orientation), management strategies (e.g., branch pruning, root pruning, fertilization, irrigation), and responses to environmental variation (e.g., climate change, soil depth and fertility, fluctuating water table, latitude). Virtual experiments could save substantial time and money by helping long-term field experiments target key interactions, place sensors in the most important locations, and focus on the most promising treatment designs. In particular, coupling Hi-sAFe with an optimization algorithm that can efficiently explore the many possible combinations of system characteristics would provide a novel, powerful tool.
While the algorithms and functionality of Hi-sAFe were primarily built around alley cropping systems, simulations of other agroforestry systems are also possible. For silvopasture systems, a perennial grass or grass-legume mix can be used as the crop and grazing events can be mimicked via repeated crop harvest and concurrent fertilization with manure. However, some complex interactions in silvopasture are not included, such as soil compaction by livestock or heterogeneous grazing based on differences in temperature, shading, or forage quality below trees. The shade and root competition effects of windbreaks or hedgerows can also be simulated by only including PBC along one axis. However, no effects of wind speed reduction, erosion control, or sedimentation and slope change are included so far. Sustainability 2019, 11, x FOR PEER REVIEW 18 of 25

Model Limitations
While deterministic and non-process-oriented models can also be an effective approach to represent complex natural systems [84,85], fully understanding agroforestry systems requires an approach that does not oversimplify their complex interactions [16,17]. Nevertheless, limitations of process-based models do exist. Results can be complex and difficult to understand [84], parameter

Model Limitations
While deterministic and non-process-oriented models can also be an effective approach to represent complex natural systems [84,85], fully understanding agroforestry systems requires an approach that does not oversimplify their complex interactions [16,17]. Nevertheless, limitations of process-based models do exist. Results can be complex and difficult to understand [84], parameter sensitivity is often ignored [86], and results can be difficult to communicate [85]. Furthermore, the development of a process-based model requires a large initial investment, and ongoing maintenance can be cumbersome [87]. Beyond its general complexity, some specific limitations of Hi-sAFe are outlined below. Overcoming some of these limitations will require new field investigations to quantify processes, while others are already understood and can be implemented now.
Hi-sAFe can only simulate deciduous broadleaf trees with crowns approximated as ellipsoids. There is also no algorithm included for nitrogen fixation by trees. Adding other crown shapes, such as paraboloids or cones, would require modifications to the light interception module and more flexible specification of aboveground allometry. Simulating shrubs or coniferous tree species would require new algorithms for phenology, carbon allocation, and aboveground allometry. Shrubs would require substantial model modifications as they have no trunk biomass pool and no DBH, which is the trait that drives the core allometric relationships.
Leveraging fruit and nut trees and shrubs within alley cropping systems has been proposed as a frontier in temperate agroforestry research [12,64]. Inclusion of a fruit compartment was initially avoided in Hi-sAFe because (1) most alley cropping experiments focus on timber trees [5] and (2) allocation of resources to fruit production is a complex, poorly understood process in many species (e.g., [88][89][90][91]). However, even if the dynamics of fruit production cannot be easily modeled, it may still be important when modeling some species to include a fruit sink that competes with structural growth for carbon and nitrogen allocation. For example, over half of annual carbon assimilation in olive trees can be allocated to fruit production [89].
While Hi-sAFe includes a wide range of tree responses to environmental stimuli, some processes have not yet been incorporated. Photosynthetic assimilation is known to be sensitive to atmospheric carbon dioxide concentration and leaf temperature [92]. STICS incorporates these sensitivities into the crop growth algorithm, but no such algorithm currently exists for trees in Hi-sAFe. Tree LUE is only sensitive to leaf age, water stress, and nitrogen stress. Furthermore, while tree phenology can show complex responses to heat, frost, and drought stress [93,94], Hi-sAFe does not include sensitivity of bud burst timing, growing season length, nor leaf fall duration to these stresses other than frost-triggered leaf senescence. Additional complexity also arises as these responses can vary with tree age in many temperate deciduous trees [95]. Phenological modeling has been largely empirical, although process-based models are being developed and could be incorporated [96,97].
Plasticity in crown architecture can play a critical role in tree-tree interactions, constituting an important mechanism of niche complementarity in multispecies systems [79]. In Hi-sAFe, crowns stop growing horizontally when they touch other crowns, but there is no anticipatory response of crown architecture to shading stimuli (e.g., increased allocation to height growth to avoid shade). Furthermore, crown interactions in Hi-sAFe are currently only possible with single-tree scenes. Inclusion of more sophisticated interactions among tree crowns will require an enhanced understanding of how trees anticipate and respond to shade, although any substantial improvements will be limited by the ellipsoid approximation of crown shape.
Belowground processes can be difficult to quantify, making process-based root and soil models hard to validate. The Hi-sAFe root growth model developed by Mulia et al. [47] is simple and flexible, though it has only been validated using potted plants. Further investigation of root architecture in field-grown trees, especially in agroforestry systems, will be necessary for a robust validation of this model. Root growth in Hi-sAFe currently occurs only while trees have active leaves, with phenology and lifespan uniform across depths. However, root growth has been observed year-round in some species, and fine root phenology and lifespan can change with depth or temperature [98].
The tree carbon assimilation algorithm is a simple empirical model, with LUE (Equation (4)) and ∆LFR* (Equation (11)) sensitive to water and nitrogen stress in a simple, multiplicative manner. Further investigation is necessary to explore alternative representations of stress responses. For example, one alternative to Equation (4) (and analogously for Equation (11)) would be to apply only the minimum of the two stresses: Furthermore, the respiration in trees is not explicitly accounted for, limiting some more complex investigations of climate on tree growth. Some limitations exist in the aboveground impacts of the trees on crop growth. Since light interception is calculated at a half-hour resolution, the impact of shorter transition phases between shade and sun on crop growth are not modeled. Recent findings suggest that this limitation could overestimate crop photosynthesis and light-use efficiency [99]. Another limitation is that while Hi-sAFe does account for the impact of tree shade on crop radiative forcing on crop temperature, the tree canopy does not modify air humidity due to the mixing of crop and tree transpiration. While this process is not major during windy days, it may reduce the transpiration demand of the crop during still days. The Shuttleworth-Wallace equations [100] could be added into the model to account for this effect, but this would substantially increase the computational requirements. This algorithm was one of the major reasons for the discarding of the HyPAR model [32]. Similarly, wind reduction by trees is currently not included in Hi-sAFe. Inclusion of a wind reduction effect would require wind direction as a new meteorological input, which is not available at many meteorological stations.
Some model limitations arise due to each cell in Hi-sAFe being managed by an independent instance of STICS. For example, while the light interception module in Hi-sAFe can handle sloping scenes, belowground water and nutrient fluxes are not responsive to slope as there are no horizontal fluxes between cells. Similarly, the decision to harvest the crop occurs independently for each cell. This can lead to very different harvest dates across cells because interactions with the tree can influence crop growth, maturation, and grain drying. In practice, however, a farmer would harvest the entire scene at one time. Incorporating this within Hi-sAFe will require enhanced communication with STICS so that harvest can be forced across all cells if, for example, a certain proportion of cells have reached the harvest stage. This limitation has less of an impact on grain crops, as harvest can occur once all areas are ready. The effect on forage crops that require multiple cuttings, however, can be large.
The STICS crop model also has some important limitations that are especially relevant to its use within Hi-sAFe. While most temperate crops have already been parameterized for the model, few tropical crops are available. STICS was also not designed to simulate the fate of carbon in deep soil layers. While over two-thirds of roots in both annual and perennial crops are typically above 50 cm [101], a substantial portion of tree root growth and decay occurs much deeper. Consequently, Hi-sAFe is currently not well suited for modeling long-term soil carbon sequestration in agroforestry systems.

Model Implementation and Distribution
Hi-sAFe version 4.0 is available online free of charge (at www1.montpellier.inra.fr/wp-inra/hisafe/en/). The model runs on PC, Mac, and Linux platforms and requires Java version 1.8 or later. A suite of tools for building, running, and analyzing Hi-sAFe simulations is also available via the hisafer R package (www.github.com/kevinwolz/hisafer). The package includes a vignette to introduce users to the available functions and an online issue tracker where users can report bugs or suggest changes.

Conclusions
In this paper, we introduce Hi-sAFe as a powerful tool for understanding complex agroforestry systems around the world. Designing a model capable of predicting long-term tree growth and crop yield, both as monocrops and as mixtures, was a significant challenge that took many years of effort. Compared to other existing agroforestry models, Hi-sAFe has several unique features, such as (1) 3D plastic tree root architecture reaction to water and nitrogen availability, (2) an integrated algorithm for water and nitrogen competition with no biased priority for trees or crops, (3) a fluctuating water table, (4) spatially explicit tree and crop management interventions, and (5) a dedicated R toolbox for designing, running, and analyzing simulations. Nevertheless, further work on Hi-sAFe is needed to parameterize tree species, update algorithms as agroforestry science progresses, enhance analytic tools, and improve computational performance. Overall, Hi-sAFe is a powerful tool that is ready to (1) extrapolate results from existing agroforestry experiments to new system designs, management strategies, and pedoclimatic environments, (2) identify and prioritize knowledge gaps, (3) generate quantitative hypotheses for testing in the field, and (4) select the most promising system designs to include in field experiments. Agroforestry models such as Hi-sAFe are critical tools for enhancing scientific understanding and, ultimately, facilitating widespread adoption of agroforestry as a sustainable land-use practice. We gratefully acknowledge all these partners for supporting our long-term efforts on this model.