Water Flow in the Soil Profile Under an Agroforestry System in the Colombian Amazon Based on HYDRUS-1D Model

1. Introduction

Soil water availability is an essential variable for hydrological understanding[1], which is influenced by factors such as environmental conditions (precipitation, evapotranspiration, radiation [2]), soil physical properties[3,4], crop management[5,6], composition and structure under agroforestry arrangements [7,8,9,10]. In this sense, the importance of soil water monitoring has been reported by different studies[11,12,13], among which the generation of management and adaptation schemes[2] in crops such as cocoa[7,14] stands out.

Regarding soil hydraulics, soil structure controls retention functions and hydraulic conductivity[15]. In this sense, knowing the spatial (i.e., at depth level) and temporal variability of water in the soil is fundamental to accurately describing the processes of infiltration and water storage[16]. Therefore, performing simulations of soil water dynamics can be a valuable tool to evaluate different land uses and their impact on water availability [17].

Developing accurate simulations requires that the hydraulic properties of the soil, i.e., hydraulic conductivity and water retention, are adequately known. These parameters can be derived from soil physical properties using pedotransfer functions from direct or inverse measurements[18], or through models[19,20], whose objective is the assessment of water availability in the soil profile to determine states of water deficit [21,22,23] that affect crops development.

Cocoa is grown under different agroforestry structures[24], which impacts the water status of cocoa trees[25,26]; therefore, knowing the dynamics of soil drainage under cocoa production systems is important for managing a proper water balance. Understanding the dynamics of soil drainage is necessary to know the water demand because when there are situations of stress caused by water deficit, physiological processes such as photosynthesis are affected[27,28,29,30], resulting in adverse effects on production[14,31,32]. According to the above and based on characteristics such as water availability in the soil profile, physical properties, hydraulic properties and precipitation, the objective of this study was to model the volumetric soil water content (θ, cm³ cm^-3) under agroforestry systems (AFS) with cocoa (Theobroma cacao) as the main crop in the Colombian Amazon. The information generated will give a better understanding of the dynamics of soil drainage in agroforestry systems in the Amazon under drought conditions, which may become more frequent and severe due to climate change[7,33].

2. Materials and Methods

2.1. Study Site and Monitoring of Environmental Variables under Agroforestry Systems

The study was conducted in a cocoa agroforestry arrangement at the Macagual Research Center - University of Amazonia-Colombia (1°37’N y 75°36’W). This region presents a warm-humid climate, characteristic of the tropical rainforest ecosystem, with an average annual rainfall of 3,800 mm, a sunshine of 1,700 hours year^-1, an average temperature of 25.5°C and relative air humidity of 84%. Four cocoa plots (20 × 50 m) were established into a cocoa agroforestry system. The agroforestry system is composed of cocoa plants planted at a distance of 3 × 3 m and timber species such as Cariniana pyriformis and Calycophyllum spruceanum with a density of 55 trees per ha^-1 that generate a transmitted radiation of 28%, exhibiting a typology classified as diversified multistrata shade[24].

Monitoring was carried out in each plot, from day of year 66 to 85 of 2015 (March 7 to 25). A WatchDog 2900ET (Spectrum Technologies, Inc., USA) meteorological station was placed under the canopy of the agroforestry arrangement at the height of 1.5 m, and the precipitation (PP mm), air relative humidity (RH_a %), air temperature (Ta, ºC), and photosynthetically active radiation (PAR, µmol m^-2 s^-1) monitored every minute. The vapor pressure deficit (VPD, kPa) was calculated from the air temperature and relative humidity recorded minute by minute following the methodology proposed by Allen et al.[34]. With the different environmental variables, the reference daily evapotranspiration (ET₀; mm day^-1) was calculated according to the FAO Penman-Monteith method[34]:

$$E{T}_0=\frac{0.408\mathsf{\Delta }\left(R_n-G\right)+\gamma \left(\raisebox{1ex}{$900$}\!\left/ \!\raisebox{-1ex}{$T$}\right.+272\right)u_2\left(e_s-e_a\right)}{\mathsf{\Delta }+\gamma \left(1+0.34u_2\right)}$$

(1)

where $R_n$ is the net radiation absorbed by the surface [MJ m^-2 day^-1], G is the ground heat flux density [MJ m^-2 day^-1], T is the air temperature [◦C], $u_2$ is the wind speed [m s^-1], $e_s$ is the saturation vapor pressure [kPa], $e_a$ is the actual vapor pressure [kPa], ∆ is the slope of the vapor pressure curve [kPa ◦ C^-1], γ is the psychrometric constant [kPa ◦ C^-1]. The climatic water balance (CWB) was calculated as PP-$E{T}_0$ according to FAO[35].

2.2. Determination of Soil Hydraulic Characteristics and Properties

The 5TE sensor (Decagon Devices Inc., Pullman, WA, USA) connected to an Em50 data logger (Decagon Devices Inc, Pullman, WA, USA) between 10 and 100 cm depth was used to monitor minute by minute the volumetric soil water content (θ, cm³ cm^-3) during the monitoring period in each plot. From the texture and bulk density at each depth evaluated and using the pedotransfer function of the ROSETTA software[20], initial appreciations were derived to estimate soil hydraulic parameters (θ_s, θ_r, α ,n and K_s Table 1). This function predicts the saturated hydraulic conductivity (K_s) and the saturation point (K₀) on the saturation curve, which may differ from the K_s match. The K₀ was used as the initial parameter value of K_s (Table 1).

2.3. Modelo Hydrus-1D

To simulate water movement in the soil profile under the agroforestry system, the program Hydrus-1D[36] was used. This program numerically solves the Richards equation for water flow[36], described as follows:

$$\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial Z}\left[k\left(\frac{\partial h}{\partial z}+1\right)\right]-S\left(Z, t\right)$$

(2)

Where: θ is the volumetric soil water content (cm³ cm^-3), t is the time (d), z is the vertical coordinate space (cm), h is the pressure height (cm), K is the hydraulic conductivity (cm d^-1), and S is the water uptake term by plant roots (cm³ cm^-3 d^-1). The hydraulic properties of the unsaturated soil were described using the van Genuchtene Mualem functional relationships[37] obtained from the ROSETTA software[20], as follows:

$$\theta \left(h\right)=\theta r+\frac{\theta s+\theta r}{{\left[1+{\left(\alpha h\right)}^n\right]}^m} h<0$$

(3)

$$\theta \left(h\right)=\theta s h\ge 0$$

(4)

$$k\left(h\right)=K_s{S}_e^l{\left[1+{\left(1-S_e^{l/m}\right)}^m\right]}^2$$

(5)

$$\mathrm{with} m=1-\frac{1}{n}$$

(6)

Where: θ_s is the saturated water content (cm³ cm^-3); θ_r is the residual water content (cm³ cm^-3); m, α and n are empirical form factors in the water retention function, where m = 1 - 1 / n; K_s is the saturated hydraulic conductivity (cm day^-1); l is the form factor in the hydraulic conductivity function and Se is the relative saturation.

2.4. Data Analysis

A descriptive analysis (means and frequencies of the variables) was carried out for the climatic variables (relative humidity, temperature, precipitation, vapor pressure deficit). Volumetric soil water content (θ, cm³ cm^-3) values obtained through simulations with HYDRUS-1D were compared with the observations made in the Theobroma cacao agroforestry arrangement, and thus the model fit was determined[38]. The assumptions of normality and homogeneity of variance were evaluated using the residuals studied. The Bayesian and Akaike information criteria (BIC and AIC) and the adjusted R² were used for model selection. For model goodness-of-fit, root-mean-square error (RMSE) and Percent Bias (PBias) were calculated, as well as modeling efficiency (EF), also known as Nash-Sutcliffe Efficiency (NSE)[25]. A standard regression was used to describe the relative relationship between the observations and the simulation (slope), which allows for identifying any lag or deviation between the simulated and observed values[39]. For this purpose, scatter plots were created showing the line of fit and the line of best fit (i.e., y=x). The model's goodness of fit was evaluated using the ggof function[40] of the HydroGOF package of R-Project version 4.1.1.[41] using the implemented interface of the R platform in InfoStat[42].

3. Results

3.1. Micrometeorological Conditions and Soil Water Content under Agroforestry Systems

During the study period, micrometeorological conditions were monitored simultaneously (Figure 1). The radiation levels ranged from 2 to 1,242 µmol m^-2 s^-1 the day with an average of 380±22 µmol m^-2 s^-1 which influenced humidity and ambient temperature. The RH_a presented minimum values of 48.1 with an average of 80.3 during the day and minimum values of 80.5 with an average of 97.2 during the night. The temperature ranged from 21.5 to 33.5 ºC, averaging 25.5±0.1 ºC). These two variables influenced the vapor pressure deficit, which ranged from values very close to zero to 2.1 kPa with averages of 0.37±0.02 kPa. Likewise, a total of 315.7 mm of water accumulated as a result of 46 rainfall events with an average rainfall of 6.86±1.75 mm in 45.9±7.81 minutes with an intensity of 0.12±0.03 mm minute^-1. There were seven moments in which ET₀ was higher than 5 mm d^-1, which coincided when the environment had higher PAR values of 1041 µmol m^-2 s^-1, very low atmospheric RH_a levels (around 45%), higher ambient temperatures of 32.2 ºC and a vapor pressure deficit of 1.6 kPa. During those moments, the soil volumetric water content levels reached values of 0.32 cm³ cm^-3 and a CWB of -5.7 mm. CWB exhibited negative values during most of the study period. The θ in the soil profile ranged from 0.32 to 0.13 cm³ cm^-3, going from 0.33 cm³ cm^-3 at 0-20 cm depth to 0.11 m³ m^-3 at 80-100 cm depth, enlightened by a linear model with a high level of fit (cm³ cm^-3 =0.387-0.005*depth (cm), P<0.0001 r²=0.86, Figure 2). The predicted values of θ in the soil profile obtained from the HYDRUS-1D model closely resembled the observed data (Figure 2).

3.2. Soil Water Properties

Significant differences were found for θ_f (P<0.0001), θ_r (P<0.02), θ_s (P<0.03), α (P<0.0001), n (P<0.0001), K_s (P<0.01), and K₀ (P<0.0018 at different depths (Table 1). Likewise, there were variations in relation to the level of compaction in the soil profile that affected water availability, finding a negative relationship with θ_s (r: -0.67 P<0.0001) contrary to the tendency between K₀ (r: 0.94 P<0.0001).

3.3. Hydrus-1D Calibration and Validation Model

The simulated data allowed explaining θ in each of the soil profiles, calculating each of the peaks that occur at the time of rainfall events, as well as the amount of θ as the depth increases (Figure 2). The relationship between observed versus predicted θ was estimated and simulated using the HYDRUS-1D model for each depth of the soil profile evaluated and the estimated values were found to be very close to the y=x line (Figure 3). Based on the goodness-of-fit estimators of the model, a small bias or deviation was found between the simulated and observed values measured from the RMSE, PBias and NSE, which were calculated at each depth of the soil profile (Figure 3). Likewise, a high coherence level was found, measured by R², which was greater than 0.96 (Figure 3).

4. Discussion

4.1. Microclimatic Conditions, water Content and Soil Properties under Agroforestry Systems

The AFS have been shown to influence environmental variables, which affect the water status of cocoa farms, as well as the distribution of water in the soil profile [25,26]. Under AFS, trees used as shade canopy regulate the transmitted radiation affecting soil evaporation and leaf transpiration[43,44] and consequently the vapor pressure deficit (VPD) and evapotranspiration (ET₀)[45,46]. These environmental drivers can lead to direct and indirect effects on plant and crop physiological functioning (stomatal conductance, gas exchange, sap flow velocity)[25,47,48], which in turn, depending on the water adjustment made by soil-plant-atmosphere, translates into a direct effect on soil water content availability[7,33,49].

Given the mechanisms and processes discussed above, our study showed that ET₀ dynamics depend mostly on soil water availability from rainfall events. Between each rainfall event, it was observed that water moved through the soil profile, and due to percolation, availability was reduced in the first centimeters of depth[50], significantly increasing percolation, which allowed a high water content to be observed in the soil profile during the entire monitoring[51]. This is due to the mechanism that relates the precipitation variability with the displacement of the location of water availability in the soil profile, which results in a disproportionate increase in the depth of water penetration with the amount of precipitation[52]. In addition, depending on the amount of water (mm) that reached the soil by precipitation, the soil presented a sensitivity in volumetric water content (θ) instantaneously, there being a short-term preferential flow during rainfall events that can generate infiltration through the soil profile[53].

On the other hand, there are soil characteristics that affect ET, which are related to bulk density (ρb) and organic carbon content[54]. For example, ρb increased proportionally to depth, however, probably due to the action of the root system ρb was reduced between 20 and 40 cm depth[55,56]. Studies have described that the cocoa crop tends to present a greater concentration of roots in the first 60 cm of depth to capture a greater amount of water[57,58], a situation that modifies the ρb. Likewise, small-scale determinations of soil texture and organic carbon content (CO%) showed different variations within the soil profile. For example, between 20 and 60 cm, a higher clay content (clay loam) was found, a condition that affected water mobility in the soil profile[13] as well as the amount of CO in the soil, which was 2% in the first 10 cm of depth, an amount that was reduced by half after 60 cm, conditions that have been described as causing the modification of water availability in the soil[59]. In addition to the above, the field capacity (θf) in the soil profile was inversely proportional to the ρb, specifically in the zone with the highest concentration of roots[50,60,61,62] as well as a relationship was also found in this way between ρb and θ_f with moisture content, as well as with clay and organic matter content that favored a higher soil saturated water content (θ_s)[63,64].

4.2. Hydrological Properties of the Soil

Soil physics (texture, bulk density) influenced the hydrological properties of the soil as well as water availability related to infiltration, percolation, surface runoff, inter (or soil) flow and groundwater flow[59,65] in the soil profile. For example, our study found that saturated water content (θ_s) presented significant differences at the soil depth level, with higher saturation at the 20-40 cm depth corresponding to the soil layer with higher root density[66]. In addition to the above, a higher θ_s was related to clay texture, which has a maximum water holding capacity and a minimum water release property[67,68]; furthermore, a significant effect of θ_s on K_s and n values was observed where small changes in θ_s promoted variations of n and K_s throughout the soil profile[69].

On the other hand, soil water retention is mainly related to changes in bulk density and its magnitude depends on soil texture[70,71]. Likewise, the soil water retention function is subject to the hysteresis phenomenon, which manifests itself as a difference between the soil wetting and drying equilibrium curves (hysteresis cycle)[72]. Based on the above, in our study, we presume that a hysteresis phenomenon is present in the soil profile due to changes in water availability in the soil profile (θ) due to multiple rainfall events; added to this, variations in soil texture (sandy clay loam and clay loam), generated a bottle effect due to the lack of uniformity in the shape and sizes of both individual pores and interconnected pore networks[73,74]. As a consequence, the parameter α, which is more related to soil structure, is more affected by the different wetting processes[75]; in addition, these conditions caused the residual water content (θ_r) to present variations in the soil profile due to higher suction (i.e., loam and clay loam soils with higher α values), as we observed in the depth of 20-40 cm[68,76], which causes almost all soil pores to remain filled with water at high tensions[72]. For the case of parameter n, it presented greater stability as it was associated with texture than with soil structure[77], presenting the lowest values in the soil with the sandiest fraction (at depths 0-20 and 60 -100 cm) and higher values in the most clayey texture (depths 20-60 cm)[68,78]. It should also be noted that the low variability of n values may result in a considerably different shape of the soil water retention curve[77].

The dynamics of soil structure and its associated factors (e.g., root biomass, root system activity, soil fauna, and organic matter input) decrease with soil depth[79,80]. This trend causes saturated hydraulic conductivity (K_s) to decrease with increasing soil depth[81]. Our study found that these factors can influence K_s in the first 20 cm depth. Moreover, the results obtained were in agreement with da Silva et al.[82] and Marques et al.[83], who, when analyzing K_s with different vegetation types in the lower Amazon, found higher values of this property at the surface, alluding that K_s is more sensitive to changes in soil properties. Also, our results indicated that a high contribution of organic matter (by the accompanying trees and cocoa) in the surface layer improves the structure and influences the hydraulic and translocation properties of clay, forming very hard layers in some depths. This effect caused that when observing the behavior of K_s by depth, a relationship between K_s decrease with depths with higher clay content was found [84,85]. In addition, we found a highly significant relationship between decreases in K_s and θ_f with increases in K₀ and ρb in the soil profile[86], as a positive correlation between Ks and θs values, which may be associated with the presence of larger pores[77].

4.3. Hydrus-1D Calibration and Validation Model

Under field conditions, soil properties exhibit variations in vertical and horizontal directions[66]. Also, soil hydraulic properties that control soil water dynamics and redistribution vary from one point to another[66]. Thus, in this study, the van Genuchten model was selected to determine soil hydraulic properties and water retention variables (θ_f, θ_s, α, n, K_s). By inverse modeling, the observed and simulated soil water content was adjusted by the Rosetta pedotransfer function[20,87]. In this sense, the results between observed and simulated data presented very similar variation trends, with values in water content (θ) estimated close to the observed data [82]. The HYDRUS-1D model was able to simulate the θ at the different depths accurately, the dynamics consistently reflect rainfall events and the data indicate that the observed soil water content is well represented by the modeling results, with R² values equal to 0.96 and 0.99[88]. However, some discrepancies were observed between the measured and simulated θ, specifically on day 81 of 2015 (March 22), indicating that the model could not effectively capture the observed data, particularly for high values. This same observation was obtained by Raki et al.[89], Silva et al.[82] and Greccoa et al.[90] where they concluded an underestimation where the model did not respond to small precipitation events that occurred in short periods, overestimating the θ values measured in the field relating them to the high infiltration rate considered by the model, as well as overestimating evapotranspiration rates, which can reduce soil water content. In addition, Rezaei et al.[91], similar to us, associated the small variations between simulated and observed values to the soil hysteresis phenomenon, as discussed above.

The simulated and observed values of soil water content presented R² that ranged from 0.96 for the 20 cm depth to 0.99 for the other depths (40, 60, 80, 100) and the RMSE values ranged in minimal values (0.00002 to 0.00003). Being the depth of 60 the one that presented the lowest RMSE (0.0002), where a higher and almost constant θ was observed[66]. In this sense, the simulated values of soil hydraulic parameters such as it θ were close to the measured ones. The above shows the capability of the HYDRUS-1D software to simulate field conditions, even under agroforestry systems [92]. Moreover, this method and the goodness of fit produced the smallest variations, thus obtaining a great prediction[93,94]. Considering the above results and the entire time series used, both simulation and validation, for use in the HYDRUS-1D model describe well the soil conditions under the cocoa-woody agroforestry system. Also, we believe that simulation based on field data can reproduce the water movement and change process more accurately[68]. Also, the impact of hysteresis on the soil water content process needs to be further considered to improve the simulation accuracy and model results[95].

5. Conclusions

When modeling the volumetric soil water content (θ, cm³ cm^-3) under agroforestry systems (SAF) with cocoa, a small bias was found as a high level of coherence between the observed values, based on the goodness of model fit measured by the Root Mean Square Error (RMSE), Nash-Sutcliffe efficiency (NSE), Percent Bias (Pbias), and coefficient of determination (R²). During the monitoring period, it was found that cocoa cultivation under agroforestry systems presented a negative climatic water balance, a water status that requires considerable attention since evapotranspiration (ET₀) conditions greater than 5 mm d^-1 were present in conditions where the radiation was greater than 1,041 µmol m^-2 s^-1, with relative humidity levels in the atmosphere of 45% and higher ambient temperatures of 32.2 ºC that generated a vapor pressure deficit of 1.6 kPa, which can be worrying in current conditions of climatic variation with a tendency to temperature increases.