Evaluating tree growth factors into species-specific functional soil maps for improved agroforestry system efficiency

Agroforestry systems play an important role in sustainable agroecosystems. However, accurately and adequately quantifying the relationships between environmental factors and tree growth in these systems are still lacking. Objectives of this study were to quantify environmental factors affecting growth of four tree species and to develop functional soil maps (FSM) for each species in an agroforestry site. The diameter at breast height, absolute growth rate (AGR), and neighborhood competition index of 259 trees from four species (northern red oak [Quercus rubra], pecan [Carya illinoinensis], cottonwood [Populus deltoides], and sycamore [Platanus occidentalis]) were determined. A total of 51 topsoil samples were collected and analyzed, and 12 terrain attributes were derived from the digital elevation model. The relationships between AGR, soil, topography, and tree size were analyzed using Spearman correlation. Based on correlation analysis, FSM for each species were generated using the k-means cluster method by overlaying correlated soil and terrain attribute maps. Results showed tree size and terrain attributes were driving factors affecting tree growth rate relative to soil properties. The spatial variations in AGR among functional units were statistically compared within tree species and the areas with larger AGR were identified by the FSM. This study demonstrated that FSM could delineate areas with different AGR for the oak, cottonwood, and sycamore trees. The AGR of pecan trees did not vary among functional units. The generated FSM may allow land managers to more precisely establish and manage agroforestry systems.


Introduction
Agroforestry systems are being increasingly adopted in North America because of their advantages in improving economic incomes, wildlife diversity, carbon storage, and mitigation of regional climate change (Jose et al. 2019). Generally, agroforestry systems consist of trees, grasses, and livestock in the same land base to achieve sustainable development goals (Chará et al. 2019). As the main component of the agroforestry system, trees, and tree growth dynamics must be monitored and understood to make the entire system more efficient and productive (Villegas et al. 2009).
The tree growth rate in any agroforestry system may be highly variable within the same species and among different species due to variability of multiple biotic and abiotic factors (Chamagne et al. 2017;Stephenson et al. 2014). Influences and interactions of biotic and abiotic factors may differ among soil type, tree species, and management. In small-scale areas where single or a limited number of tree species are planted, micro-topography variation has been found to strongly affect soil nutrient distribution, water distribution, and biological processes, which may influence variations in individual tree growth rates (Bassett 1964;Cavelier et al. 2000;Kübler et al. 2020;Rhoades 1996). Additionally, the individual tree growth rate may vary across different sizes and ages of trees (Dold et al. 2019;King et al. 2006). Therefore, quantification of relationships among tree genetics, tree growth, soil variability, and soil-landscape features are needed to maximize efficiency in agroforestry systems.
The development of pedometrics and functional soil mapping techniques provides an opportunity to quantitatively explore soil-landscape relationships for a specific agricultural purpose (Owens et al. 2020;Jiang et al. 2021). In functional soil maps (FSM), areas with similar soil, terrain, and landscape attributes can be classified as functionally homogeneous areas based on fuzzy logic classification methods (Ashtekar and Owens 2013;Owens et al. 2020). Generated soil maps can then help land managers make better land use and management decisions for different functional areas to improve land productivity and achieve sustainable development simultaneously. Previous studies have primarily focused on exploring uses of FSM in precision agriculture and land management (Salahat et al. 2012;van Zijl et al. 2014). Penney et al. (1996) developed FSM for the improvement of fertilizer management for spring wheat (Triticum aestivum L.). Bakhsh et al. (2000) investigated the relationship between soil attributes and corn (Zea mays) -soybean (Glycine G. max) yield variability to delineate areas with different yields using FSM. Zhu et al. (2013) generated FSM to capture the variations of soil moisture and crop yield to support precision farming techniques. However, the study of FSM in agroforestry systems and the applications on precision tree management have rarely been studied to date. Agroforestry systems can benefit from FSM techniques that can quantitatively analyze effects of environmental factors on tree growth within agroecosystems. In this study, FSM techniques integrated with soil and topographical attributes are presented to capture the spatial variations of tree growth and generate functional soil property maps for four tree species. The main objectives of this study were to (i) investigate the variation of growth rate for four species of trees in an agroforestry system, (ii) explore the effects of soil, topography, and tree size on the spatial variability of growth rate for four species of trees, and (iii) develop and evaluate the species-and site-specific FSM for explaining the spatial variation of tree growth rates. It was hypothesized that soil and terrain attributes could influence the tree growth and FSM could capture the spatial variation of tree growth rates for use in precision land and tree management.

Study area and trees
The study site was located in the University of Arkansas Agricultural Research and Extension Center in Fayetteville, Arkansas (latitude 36.1°, longitude -94.17°). It is a 4-hectare area where the main landform is within a sub-basin of the Illinois River watershed in Ozark Highlands. The elevation of the study area ranges from 378 to 386 m a.s.l. with slope varying from 0 to 9° (Fig. 1). The mean annual temperature of the study area is 13.7°C and the mean annual precipitation is 1217 mm (Sauer et al. 2015). Soils of the study area were derived from colluvium and residuum, which can be classified into fine-silty, siliceous, active, mesic Typic Fragiudults; fine-silty mixed, semiactive, thermic Typic Paleudults; loamyskeletal, siliceous, active, mesic Glossic Fragiudults; and fine-silty, mixed, active, mesic, Aquic Fragiudults (Harper et al. 1969).
The study site was established from ungrazed pasture in 1999. The trees consisted of northern red oak (Quercus rubra L.), pecan (Carya illinoinensis K. Koch), cottonwood (Populus deltoides W. Bartram ex Marshall), and sycamore (Platanus occidentalis L.). Trees were planted in 1999 (pecan), 2000 (oak), and 2014 (cottonwood, sycamore). The distance between tree rows was 15 m. The distance between individual trees within the same rows during this study (all species except pecan had been thinned) was on average 7 m for oak (due to mortality and thinning) and 9.1, 4.6, 4.6 m, respectively, for pecan, cottonwood, and sycamore. The location of individual trees in the study area was recorded by hand-held global positioning system (GPS) receivers with 0.5 m accuracy. Forage crops of orchardgrass (Dactylis glomerata L. 'Tekapo') and big bluestem (Andropogon gerardii L.) were seeded in alleys between tree rows (15 m wide) in 2015.

Tree growth measurement
The diameter at breast height (i.e., 1.37 m above soil level) (DBH) of the trees in the study area were measured annually to represent the tree growth. The DBH was measured using steel dendrometer bands (Thomas et al. 2020). The DBH of pecan trees was measured annually from 2004 to 2019, oak from 2005 to 2019, and of cottonwood and sycamore trees from 2016 to 2019. The dead, thinned, or replaced trees during the study periods were excluded from the analyses in this study. A total of 81 pecan trees, 71 oak Fig. 1 The study area and distribution of trees trees, 50 cottonwood trees, and 57 sycamore trees remained for the analysis in this study (Table 1). The annual increment in DBH of each tree was used to calculate the absolute tree growth rate (AGR) to compare the growth variances of trees within each species [Eq. (1)] (Chi et al. 2015).
where AGR represents the absolute tree growth rate, D1 and D2 represent the first and last DBH measurements (cm), respectively. T2 and T1 represent the years of the first and last DBH measurements.
Environmental factors related to tree growth To explore the impacts of environmental factors on tree growth in the study area, terrain attributes, soil properties, and tree properties (i.e., initial tree DBH and neighborhood competition index) that related to the analysis of tree performance were selected based on previous studies (Brown and Greenberg 1998). Redundant variables from the pool of environmental variables that had large variance inflation factor (VIF) (i.e., VIF [ 10) were removed, and a total of 14 terrain and soil variables were retained in the subsequent analysis.

Soil properties
Soil samples were collected using a stratified sampling technique based on terrain attributes and were taken in alleys between tree rows (grass alleys). The number of soil samples were derived to provide multiple samples within each combination of terrain attribute, forage and tree species following previous sampling schema Gurmessa et al. 2021). The georeferenced sampling locations covered all topographic positions and were located by a hand-held GPS receiver with 0.5 m accuracy. A total of 51 topsoil (0-15 cm) samples (i.e., samples were collected using a 3.3 cm diameter core) were collected and transported to the lab. Soil samples were air-dried for 7 days at approximately 20°C, ground, and sieved through 2 mm before laboratory analyses. Sand and clay content of the fine-earth fraction were measured using the pipette method (Gee and Or 2002). Soil bulk density was measured from the core samples after oven-dried at 105°C for 72 h (Blake and Hartge 1986). Soil organic carbon content was measured by the VarioMax CN analyzer (Elementar Americas, Mt Laurel, NJ) using the high temperature combustion method (Nelson and Sommers 1996). Soil potassium, calcium, and aluminum content were measured by the inductively coupled argon plasma spectrometry (ICP, Agilent Technologies, Santa Clara, CA) from the 1:10 soil-extractant solution (Houk et al. 1980).

Digital soil mapping
Soil properties were mapped across space following the similarity model proposed by Zhu (1997). The similarity model applies rules, optimality curves, and fuzzy logic to measure the degree of similitude (i.e., membership) between the soil-forming conditions at a given location and soil-forming conditions represented by a set of prescribed soil taxonomic classes.
The similarity model has been integrated into the SoLIM (Soil-Land Inference Model) software. The soil properties and their maps were derived from the similarity model using the SoLIM software. *iDBH represents the initial diameter at breast height Functional soil maps for each tree species were generated using the k-means clustering algorithm by overlaying and grouping the soil and topography variables related to the tree growth in SAGA GIS (Hartigan 1975). The K-means clustering method is one of the most widely used unsupervised classification algorithms, which can subdivide the sampling points into several clusters with optimized clustering criteria. The sampling points within the same cluster have the smallest variances.

Neighborhood competition index (NCI)
Neighborhood competition has been shown to influence tree growth, especially for the same species planted in the same area where the trees have similar resources for growth (Chi et al. 2015). Neighborhood trees within the circular influence zone of the focal trees could be assumed to be competitors. In this study, the Hegyi index for each tree was calculated to represent the neighborhood competition index (NCI) within a defined circular zone (García 2014). The NCI of each tree within the same row was calculated using the Siplab package in R 3.6.2 [Eq. (2)].
where D i represents the DBH of the focal tree i, D j represents the DBH of competitor trees, R ij represents the distance between trees i and j.

Data analysis
All of the statistical analyses in this study were conducted using the R 3.6.2. The 'stats' and 'pspearman' packages in R were used for t-test and Spearman analysis, respectively (Savicky 2009). The independent sample t-test was used to examine whether significant differences (p \ 0.05) in relative tree growth exist within and among the four species and different map units for the four tree species, respectively. The Spearman rank correlation analysis was performed to explore the relationship between terrain attributes, soil properties, tree properties, and tree growth in the study area with a significance alpha level of 0.05 (p \ 0.05).

Variations of terrain attributes
The spatial variations of terrain attributes are shown in Fig. 2. Generally, the study area could be classified into the convergence flat bottom area (i.e., the channel area, and the south flat area), the middle steep slope area (i.e., the south area), and the upper flat area (i.e., the north and northeastern areas) per the spatial distribution of indices of slope, MRRTF, and MRVBF. The higher indices of SAGAWI, MRVBF, and FlowAcc in the southern and the channel areas indicated that these areas are wetter than the northern and northeastern areas (Burt and Butcher, 1985). A depression channel from the north to south in the middle of the study area was also observed according to the spatial distribution of FlowAcc and MRVBF (Adhikari et al. 2018). The flow accumulation in the channel area ranged from 1000 to 21,537, which was much higher than the other areas (Fig. 2).

Variation of soil properties
The spatial variations of soil properties are shown in Fig. 3. In general, the southern area (i.e., pecan area) of the site had the highest contents of soil moisture, sand, and potassium, which were 16, 22, and 27% higher than those in the northern area (i.e., oak area). The southern area had the relative lower elevation and higher SAGAWI, MRVBF, and FlowAcc that might accumulate water from the upper slope area (Zhu et al. 2013;Ashworth et al. 2021). The northern area had the highest soil organic carbon (SOC) content, which were 17% higher than in the southern area. The northern area was used as graze pasture that might cause the higher SOC content in this area Gurmessa et al. 2021). The northwestern area had the highest contents of silt, clay, calcium, aluminum. Adhikari et al. (2018) reported that the distribution of the soil properties was influenced by the terrain attributes (i.e., SAGAWI, MRVBF, FlowAcc, and slope etc.) within the studied agroforestry site. The areas with higher indices of SAGAWI, MRVBF, and FlowAcc had a higher probability of converging soil properties.
Comparisons of tree growth rate and neighborhood competition index among species The growth of the individual trees presented large spatial variations within the same species (Fig. 4). Pecan trees had the lowest averaged AGR among the four species, while the cottonwood trees had the highest (Fig. 4). Previous study had been reported that the cottonwood had higher intrinsic growth rates than pecan trees within similar condition (Chi et al. 2015). Besides, the species of earlier stages of succession have greater growth rate (Dold et al. 2019). On average, the annual DBH growth rate of oak trees (n = 71), pecan (n = 81), cottonwood (n = 50), and sycamore (n = 57) were 1.9, 1.5, 3.9, 2.3 cm, respectively. It has been reported that young cottonwood trees had a higher growth rate than other species planted in the study area (Dold et al. 2019). The growth of individual trees presented large differences within the species in the study site (Fig. 4). It appeared that the largest differences in AGR occurred in the cottonwood trees (SD = 0.91), followed by the oak trees and sycamore trees. The pecan trees had the smallest differences in AGR (SD = 0.21). The NCI of the individual trees also varied among and within species. The cottonwood and sycamore trees had a larger average NCI value than oak and pecan trees, which could relate to the neighborhood distance between the trees which was 2.3 m at planting and 4.6 m after thinning. Besides, the individual cottonwood and sycamore trees also had larger NCI differences within the species that could be related to the larger differences in DBH of induvial trees (Table 1; Fig. 4).

Effects of topography, soil, and tree size on tree growth
The influences of terrain attributes, soil properties, and tree size and neighborhood competition on the variations of tree growth rate in the study area were evaluated. Generally, tree growth rates were more Fig. 2 The spatial distribution maps of the terrain attributes in the study site. Slope slope percent, SAGAWI system for automated geoscientific analysis wetness index, MRRTF multi-resolution ridge top flatness, MRVBF multiresolution index of valley bottom flatness, VDChn vertical distance to channel network correlated with terrain attributes and tree size than soil properties based on the Spearman correlation coefficient (r) (Fig. 5). Terrain attributes could drive the soil moisture and soil nutrient variation and subsurface flow in the steeply sloping areas (Burt and Butcher 1985;DeFauw et al. 2014), and tree growth is Fig. 3 The spatial distribution maps of the soil properties related to tree growth in the study site. The soil properties maps were derived from the laboratory analysis data of the soil samples and mapped using the ordinary kriging method. SOC soil organic carbon Fig. 4 The annual absolute growth rate (AGR) of diameter at breast height (a) and neighborhood competition index (NCI) (b) for the four species trees in the study area. The upper and lower outliers represent 5 and 95% of the mean, the horizontal lines within the boxes represent the median values, the dots represent the AGR and NCI of trees, the lowercase n represents the number of the trees for each species. SD represents the standard deviation. Lowercase letters (a, b, c, and d) represent the significant differences of AGR among the species interrelated with soil moisture dynamics (Choi et al. 2007;Mathys et al. 2014). Even with micro-landscape variation, the terrain attributes could be the source of the dominant influence on tree growth (McNab 1989). Therefore, the influences of terrain attributes on the individual tree growth rate might be more prominent than soil properties in the study site.
The impacts of terrain attributes on the individual trees' growth rate also varied among the different species. For the oak trees, tree growth rates significantly increased with the indices of VDChn and decreased with the increase of indices of SAGAWI and MRVBF. On the contrary, the cottonwood tree growth rate significantly increased with the indices of SAGAWI and MRVBF. The growth rate of pecan trees was significantly related to the SAGAWI and MRRTF. The pecan and cottonwood areas (i.e., lower slope area) had greater potential soil moisture compared to other areas. For the sycamore trees, only sand content had correlations with the growth rate. The stronger effects (r [ 0.41, p \ 0.05) of terrain attributes on the individual tree growth of cottonwood were observed. Cottonwood trees grew faster in the area with greater indices of SAGAWI and MRVBF, which indicated that the area with higher soil moisture was more favorable for cottonwood growth. Conversely, oak trees in the study site grew faster in areas with a lower SAGAWI and MRVBF. Oak trees grow in various soils, but the best growth conditions are on deep and well-drained soils (Cogliastro et al. 1997).
The initial diameter at breast height (iDBH) was positively correlated (r [ 0.38, p \ 0.05) to tree growth rate (Fig. 5), and NCI was negatively correlated (r \ -0.34, p \ 0.05) with tree growth rate for the four tree species, which indicated that tree growth rate increased with the iDBH and decreased with the NCI for all of the trees in the study site. Similarly, Chi et al. (2015) reported that trees with greater iDBH grew faster than the smaller trees. These results indicate that the ontogenetic and neighborhood effects also play an important role in the tree growth at the study site. The larger trees had a higher probability to receive the illumination in the crown parts, which could be more competitive than smaller trees (King et al. 2006).

Functional soil map for different tree species
The environmental variables that significantly correlated with tree growth rate were considered in generating FSM for each tree species to capture the spatial variation of tree growth rate and facilitate the management of the agroforestry system. The total planted area for each tree species was divided into four functional map units based on the terrain attributes and soil properties related to the tree growth using the clustering method (Fig. 6). The tree growth rates in different map units from the FSM were summarized and compared in Table 2. The functional map delineated as A1 was the most favorable growing area within the oak area, the B3 and B4 units were the most favorable growing areas for cottonwood, the C1 and C2 units were most favorable for sycamore trees, and the D2 unit was most favorable for pecan trees (Fig. 6). Results of this study indicate that the functional map could accurately capture the variations in tree growth rate for the four tree species in the study Fig. 5 Spearman correlation coefficients between absolute tree growth rate and terrain attributes, soil properties, and tree properties for the four species trees in the study site. Slope slope percent, FlowAcc flow accumulation, SAGAWI system for automated geoscientific analysis wetness index, MRRTF multiresolution ridge top flatness, MRVBF multiresolution index of valley bottom flatness, VDChn vertical distance to channel network, SMois soil moisture content, BD soil bulk density, Sand soil sand content, Clay soil clay content, SOC soil organic carbon content, K soil potassium content, Ca/Al the ratio of soil calcium and aluminum contents, iDBH initial diameter at breast height, NCI neighborhood competition index site, and distinctly delineate the corresponding units with different growth rates.
The tree growth rate in A1 unit was greater (p \ 0.05) than the oak trees in other units. The optimal environment for oak growth was deep, well drained soils, located at the middle or lower slope (Sander 1990). The A1 unit had the lowest indices of SAGAWI, MRVBF, and slope, where soils would potentially be well-drained and deeper than the other oak units (Böhner and Selige 2006). Moreover, the fewer variations of terrain attributes for the A1 unit also revealed a relatively stable environment that was Fig. 6 The functional soil map (four functional map units) of the four species trees planted area in the study site. a Oak planted area, b cottonwood planted area, c sycamore planted area, and d pecan planted area. The functional map was derived from the terrain attributes and soil properties that significantly related (p \ 0.05) to the tree growth rate in the study area more favorable for oak growth. The highest cottonwood tree growth rate was observed in the B3 and B4 units. The larger indices of SAGAWI and MRVBF and lower slope for B3 and B4 unit areas indicated higher soil nutrients and moisture retention in the convergence area. Moreover, cottonwood grew best within the streamside areas with lower slopes (Brown and Greenberg 1998). Therefore, the growth of cottonwood in the B3 and B4 units was faster than the neighboring divergence areas with higher indices of slope and MRRTF. Sycamore trees in the C1 and C2 units had a significantly greater growth rate than the trees in the C3 unit. Soil texture is very important to sycamore growth (Jensen et al. 2008). Sycamore trees grew faster in the well-drained soils with sufficient nutrients and water availability (Weber-Blaschke et al. 2008). The faster growth rate could partly explain the reason soils with lower sand content had lower sycamore growth rates in the study site.
There were no significant differences in growth rate for the pecan trees from the different map units in the study site. The pecan trees were planted in the southern area which had the highest indices of SAGAWI and MRVBF compared to the areas. Additionally, the lowest slope also indicated a flattened topography in this area (Böhner and Selige 2006). Therefore, soil nutrients and moisture presented relative lower variations that might explain why the growth rates of individual pecan trees did not present large variations.
The development of FSM for site-specific purposes (i.e., tree growth management in this study) has been demonstrated to have potential for promoting fertilizer management, interpreting spatial variation of soil properties, and delineating different management zones for resource-efficient and precision management (van Zijl et al. 2014;Zhu et al. 2013). Several studies have shown that FSM could partition the field according to the soil functional characteristics and topography attributes, and capture spatial variation of soil moisture and crop yields in an agricultural landscape (Adhikari et al. 2018;Bakhsh et al. 2000;Zhu et al. 2013;Kharel et al. 2021). In this study, each functional map unit represents the area that had a similar tree growth rate and environment variables compared to neighboring units. Results also showed that FSM can effectively delineate areas with different tree growth rates based on related environmental factors (Fig. 6). Additionally, results of this study could provide site-specific maps for precise management of resources in agroforestry systems. Agroforestry producers could develop the specific tree management strategies for each species according to FSM to improve the tree growth rate. However, this

Conclusion
This study investigated growth rate variations of four tree species in an agroforestry system. Cottonwood trees had the largest and pecan trees had the lowest variations of growth rate. The impacts of topography, soil properties, and tree properties on the growth of the four species of trees were compared and analyzed. Terrain is a surrogate for the redistribution of water, which is likely a driver for differences observed in tree growth response variables, which can aid in agroforestry site management. In this agroforestry study site, soil properties, terrain attributes, initial tree properties, and tree species influenced tree growth rate variations. The terrain attributes had more prominent impacts on tree growth rate variations than soil properties at the scale of each tree species planted area. Additionally, the influence of terrain attributes on tree growth rate variation varied within species. In this study, the terrain attributes and soil properties that significantly related to the four tree species were used to generate FSMs for areas planted to oak, cottonwood, sycamore, and pecan trees. The functional units with well-drained soils (A1, C1, and C2) were mapped as the fastest growth areas for oak and sycamore. The functional units located in convergence area (B3 and B4) were mapped as the fastest growth areas for cottonwood. Results indicated that the FSMs developed in this study were able to identify and classify functional units with different tree growth rates, which is considered a promising site-specific tool for tree management in agroforestry systems.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.