Prediction model for sap flow in cacao trees under different radiation intensities in the western Colombian Amazon

In this study, we measured diurnal patterns of sap flow (Vs) in cacao trees growing in three types of agroforestry systems (AFs) that differ in the incident solar radiation they receive. We modeled the relationship of Vs with several microclimatic characteristics of the AFs using mixed linear models. We characterized microclimatic variables that may have an effect on diurnal patterns of sap flow: air relative humidity, air temperature, photosynthetically active radiation and vapor pressure deficit. Overall, our model predicted the differences between cacao Vs in the three different AFs, with cacao plants with dense Musaceae plantation and high mean diurnal incident radiation (HPAR) displaying the highest differences compared to the other agroforestry arrangements. The model was also able to predict situations such as nocturnal transpiration in HPAR and inverse nocturnal sap flows indicative of hydraulic redistribution in the other AFs receiving less incident radiation. Overall, the model we present here can be a useful and cost-effective tool for predicting transpiration and water use in cacao trees, as well as for managing cacao agroforestry systems in the Amazon rainforest.


Scientific Reports
| (2021) 11:10512 | https://doi.org/10.1038/s41598-021-89876-z www.nature.com/scientificreports/ pressure deficit (VPD) and photosynthetically active radiation (PAR) are the microclimatic variables that have the largest impact on sap flow; however, this assumption may be an over-simplification 16,17,34 . In a cacao AF, the shade of companion species impacts many other microclimatic variables such as incident radiation, diurnal amplitude and mean temperature, and relative air humidity. These microclimatic variables may also impact the diurnal production and transpiration of the crop, having effects on the crop that could be independent of the ones caused by VPD or PAR alone. In such cases, accurate sap flow predictions require more complex models that take into account the simultaneous effects of these key microclimatic variables, which are affected by the modification of the AF structure. Once these models are developed, they can be used to determine diurnal patterns of tree sap flow in each type of AF; posterior calculations of sap flow would then only require the measurement of these microclimatic variables. These models would ultimately be more affordable and practical tools to calculate transpiration and water use of single cacao trees [15][16][17][18]33 as well as of the entire cacao plantation grown in different AFs. Blooming in cacao trees occurs during the relative dry season, and flower and fruit set productions would depend on plant water status and C assimilation rates 35 . Diurnal sap flow rates are a proxy of diurnal tree transpiration rates, gas exchange and C fixation 36,37 that, during the onset of the reproductive period, can be assigned to flower and fruit production. Thus, those AFs displaying adequate cacao sap flow dynamics during cacao reproductive period would increase the production of flowers 38 and fruit set 39 .
The overall aim of the study is to propose a cost-effective model to predict sap flow in cacao trees based on multiple microclimatic variables relevant to cacao AFs in the continental Amazon region during the relative dry period when cacao trees flower and fructificate. The climate and overall environmental conditions of cacao plantations are similar in this region, and all plantations use similar companion species to arrange the agroforestry structure. More specifically, in this study we measured diurnal cacao sap flow (V s ) and modelled it based on different microclimatic variables (air relative humidity, air temperature, photosynthetically active radiation, and VPD). We then analyzed the diurnal patterns of sap flow movement and day-night directions of V s in three cacao AFs in the Continental Amazon region. These AFs differed in the density and type of companion species, generating different diurnal incident radiation intensities among AFs. After the validation of the Vs model, we discussed the interest of this modeling approach, specifically in relation to its usefulness for calculating diurnal water use of cacao plantations under different AF managements, as well as its potential effects on cacao production.

Materials and methods
Study site and AF characteristics. Sap flow measurements were taken from cacao trees in three agroforestry systems at the Centro de Investigaciones Amazónicas (Amazonian Investigations Center) CIMAZ Macagual-Universidad de la Amazonia (137′ N and 7536′ W at 360 m a.s.l.), Colombia. The three AFs differed in structure (number and type of companion species) and incident radiation (see below). The climate is warmhumid, characteristic of Amazonian tropical rainforest ecosystems, with a mean annual precipitation of 3443 mm, 1200 sunshine hours per year, a temperature between 28.5 and 31 °C and a relative humidity between 81 and 88% (Fig. 1). Sap flow was measured over the course of two weeks during the period of minimum rainfall in the region (November to February). We selected this period for measurements because cacao plant production, transpiration, flowering and fruit-set differed the most across AFs during this relatively dry season, and Climate diagram of the Continental Amazon region, including monthly precipitation, mean monthly temperature, and air relative humidity (values calculated based on the last 30 years of climatic data). The driest month is January (high air temperature, low precipitation, and low air relative humidity). www.nature.com/scientificreports/ adequate cacao response to potential drought spells is of great economic importance. Anyhow, as major abiotic characteristics (climate, soil and topography) were equal in all AFs, and the period studied was a relatively mild dry period, we did not expect major differences in soil water contents across AFs, and expect that major differences in cacao physiological responses across AFs would be caused by the different microclimatic conditions and incident radiation among AFs. During the rest of the year, the climate is particularly humid and cloudy, and a preliminary analysis of cacao tree physiological responses across AFs did not render significant results. The three AF plots analyzed in this study were 25 × 50 m in size. The cacao plants were all clones (and thus of the same age), planted in the AFs in October 2012 in a regular pattern with a distance of 3 m between plants, irrespective of the AF. Although the cacao plantation was identical in the three plots, each AF differentiated itself based on the type and density of the companion species (henceforth called vegetation), which were planted in 2008 in rows with a north-south orientation and created an upper canopy with varying levels of shade. Two of the AFs included companion timber species (Cariniana pyriformis, Calycophyllum spruceanum), and a third AF included Musaceae species (plantain; Musa paradisiaca). The three cacao AFs compared in this study were ( Fig. 2 42 . Each sensor consists of two thermocouples spaced 15 mm apart and inserted in the main stem at 5 mm and 20 mm depth behind the cambium. In order to calibrate using a reference value of the sap flow velocity (zero), we used the procedure suggested by Burgess et al. 42 which consists of cutting the stem below the measuring needles and recorded sap velocity (V h ) afterwards. All corrections for wounds and misalignment of the probes were made according to Burgess et al. 42 . The velocity of the heat pulse (V h ) can be calculated as a function of different parameters as: where: V h is the heat pulse velocity (m s −1 ), k is thermal diffusivity of green (fresh) wood (with a default value of 0.0025 cm 2 s −143 ), x is the distance between the heater and each of the thermocouples needles (cm; 0.6 cm for this particular sensors), and T 1 and T 2 are the temperatures (°C) in the position of the needles, located at equidistant points downstream and upstream from the heater. Heat pulse velocity (V h ) was converted into sap flux density (V s sap flow rate L h −1 ) using the equation proposed by Burgess et al. 42 , as: where b is the basic density of the dry sapwood (i.e. dry sapwood weight divided by its green volume), s is the sap density (assumed to be the density of water), m c is the water content of the fresh sapwood, C w is the specific heat capacity of the dry wood matrix (1200 J kg −1 K −1 at 20 °C 44 ), C s is the specific heat capacity of the sap (assumed to be that of water, 4186 J kg −1 K −1 at 20 °C 44 ), and S is the cross-sectional area of the conducting sapwood. Raw heat pulse velocity data collected by the SFM1 Sap Flow Meters were imported into the Sap Flow Tool software (ICT International, Armidale, NSW, Australia) that calculated sap flow, diurnal sap flow and cumulative sap flow rates (L h −1 ) using the above mentioned functions.

Monitoring microclimatic variables in the AFs under study. A WatchDog 2900ET weather station
(Spectrum Technologies, Inc., USA) was placed under the canopy at a height of three meters, and moved within the plot every day to monitor the microclimatic parameters: air relative humidity (RH a ), air temperature (T a ), and photosynthetically active radiation (PAR) at a frequency of once per minute. The vapor pressure deficit (VPD) was calculated based on air temperature and air relative humidity, which were recorded minute by minute using the methodology proposed by Allen et al. 45 which records the maximum and minimum temperature and humidity values during a given period.
Modelling sap flow under different radiation intensities. The model. We first compared the microclimatic variables across AFs using linear mixed models (LMM) where RH a , T a , PAR, and VPD were the re- www.nature.com/scientificreports/ sponses variables, AFs was the treatment factor, and day and plant effects were considered random factors. Mean values were compared using Fisher's LSD post-hoc test (α = 0.01). The assumptions of normality and homogeneity of variance were evaluated using an exploratory analysis of residuals. Sap flow models were then adjusted using the totality of the data recorded during the monitoring period (n = 3816 measurements of sap flow in the four cacao plants of each AF). To predict V s (sap flow, dependent variable) based on microclimatic variables (predictor variables), new LMM were adjusted. Type of AF was included in the model as a dummy variable, and RH a , T a , PAR, and VPD as predictors. The model for each AFs was: Bayesian Information Criterion (BIC) and Akaike Information Criterion (AIC) were used as criteria for the selection of the best model. Partial residuals were studied to verify the type of relationship (i.e., linear or quadratic) to be included between each microclimatic variable and the response variable. The adjustments to the LMM were made using the function lme 46 in the package nlme of R-Project version 3.6.0 47 and using the implemented interface of the R platform in InfoStat 48 .
Validation. A simple linear regression analysis was conducted between predicted values, generated using the best (AIC/BIC) model for sap flow, and a validation dataset (n = 1740). The R 2 value was considered a measure of the predictive value of the model. We used the methodology of Vezy et al. 49 to evaluate the goodness of fit of the model. In order to evaluate the mean precision of the model with the same units as the variable of interest, we calculated RMSE and bias. The modeling efficiency (EF), also known as Nash-Sutcliffe efficiency (NSE), was also used to describe how well the graph of observed versus simulated data fits in the identity function. A standard regression was used to describe the relative relationship between the observations and the simulation (slope), which allows us to identify any lag or deviation between the simulated and observed values. For this purpose, scatter plots were created showing the fit line and the best fit line (i.e., y = x). The goodness of fit of the model was evaluated using the function ggof 50 in the package HydroGOF of R-Project version 3.6.0 47 and using the implemented interface of the R platform in InfoStat 48 .

Sap flow variability under different radiation intensities.
In the H PAR agroforestry system, the cacao plant exhibited increasing sap flows throughout the entire monitoring period, even at night (Fig. 3a). The maximum sap flow average values were 0.27 ± 0.03 L h −1 during daylight hours and 0.0300 ± 0.0023 L h −1 during the night. The AFs M PAR and L PAR , however, exhibited negative sap flows at night, with average values of − 0.0047 ± 0.0085 L h −1 and − 0.0314 ± 0.004 L h −1 , respectively (Fig. 3a). For M PAR and L PAR , nocturnal averages were significantly different from zero (t-test; p < 0.01), which may indicate a differential nocturnal hydraulic redistribution across these agroforestry arrangements. When comparing sap flow under different intensities of sun radiation, sap flow was higher in H PAR compared to the treatments, with maximum values recorded at noon. Predicted values obtained from our LMM model closely resembled the observed data (Fig. 3a,b; see below).
Monitoring microclimatic variables in the AFs under study. The air relative humidity was, on average, 90% with a variation range between 49 and 100% ( Table 2). The mean air temperature was 21.4 °C at night and 28.1 °C during the day. The mean PAR value during the day was 708.9 µmol m −2 s −1 , with a maximum value of 2310 µmol m −2 s −1 . Mean VPD was 0.28 kPa, with a maximum value of 1.86 kPa. As illustrated in Table 1, the greatest differences between AF arrangements in all microclimatic variables occurred at 13:00 h (solar time) and were significant in every case (p < 0.01); differences between arrangements in all microclimatic variables were also recorded at 10:00 and 16:00 h (p < 0.01). In the case of the L PAR arrangement, the average temperature at 13:00 h was 26.6 °C-the minimum value recorded in any of the three intensities. On average, temperatures in M PAR and H PAR were 1.81 °C and 3.12 °C higher, respectively (Table 1). There were also differences in all microclimatic variables during the night, and these differences were significant in all cases (p < 0.01); specifically, at 04:00 h, there were significant differences for T a and VPD between agroforestry systems.
For the remaining variables, differences between agroforestry arrangements followed similar patterns. The lowest values of transmitted PAR and VPD were recorded in L PAR , where there is a denser companion tree canopy; the highest values were recorded in H PAR , and intermediate values in M PAR . Relative humidity followed the opposite pattern (maximum in L PAR , minimum in H PAR ). We noticed that the PAR values at 13:00 h in the H PAR arrangement (2279 µmol m −2 s −1 ) had an mean value higher than those values obtained in the other two arrangements (Table 2; Fig. 4). When relating the sap flow values measured in the different plots to the different microclimatic variables, we found a negative correlation between V s and RH a , and a positive correlation between V s on the one hand, and PAR, T a and VPD on the other (all had correlation coefficients greater than 0.73 in general terms, and greater than 0.78 within each plot, ranked r(V s ,PAR) > r(V s ,VPD) > r(V s ,T a )) ( Table 3).

Modelling sap flow (V s ) under different radiation intensities.
The general model which estimated the relationship between observed versus predicted V s for all three AFs, showed that the estimate was very close to the y = x line (Fig. 5a-c). The model allowed us to explain the diurnal behavior of sap flow, which was affected by the microclimatic variables within each agroforestry arrangement that, in turn, impacted the incident radiation intensity that reached the plot. Specifically, the model predicted the water use of the cacao crop as a function of time, replicating phenomena such as hydraulic redistribution (HR). Indeed, V s was predicted well by the statistical model (R 2 = 0.98) that included the microclimatic regressors and exhibited significant differences between the three different AFs (Tables 4, 5). Significant differences were found between the slopes of all  , T a , PAR, VPD), with values of p < 0.01 except for T a between H PAR and M PAR (p > 0.05). Differences in intercepts were also observed depending on the AFs (p < 0.01). When analyzing the behavior of V s and the process of HR, we found that the microclimatic variables that determine this phenomenon were the RH a and VPD. This was due to the negative slopes in the model for both RH a and VPD in the radiation intensities M PAR and L PAR , which were statistically different from H PAR . An analysis of the predicted values versus the observed values produced a coefficient of determination of R 2 = 0.95 for H PAR (Fig. 5a), 0.96 for M PAR (Fig. 5b), and 0.98 for L PAR (Fig. 5c). All R 2 values were highly significant (p < 0.0001). These results indicate that the proposed modelling approach predicted the behavior of sap flow fairly well in relation to microclimatic variables of different radiation intensities (i.e., agroforestry arrangements), with the coefficient of determination being equal or greater than 0.95 for all AFs (Fig. 5). In all AFs (or radiation intensities), the regression line was a bit below the line y = x for high V s values and a bit above it for low V s values (Fig. 5). When analyzing the relationship between observed and fitted values for the model, the RMSE ranged from 0.01 to 0.02 L h −1 for the different radiation intensities or AFs (Fig. 5). Likewise, the slope of the regression of the predicted vs. observed values in the different radiation intensities was statistically different from the reference slope (y = x) for all AFs (Fig. 5). This means that, in general terms, the model overestimates the values of V s by an average of 2.8%.

Discussion
This study presents a model that empirically determines the relationship between different microclimatic variables (PAR, T a , RH a , and VPD) and diurnal sap flow patterns in cacao plants under different agroforestry arrangements in the continental Amazon region. These arrangements mainly differed in their intensity of incident radiation on the cacao crop. Incident radiation had a direct impact on microclimatic conditions that greatly differed between AFs. There are few models that include a set of microclimatic variables to predict diurnal patterns of cacao sap flow. Studies such as those developed by Köhler et al. 17 in Indonesia have also found that different microclimatic variables have significant effects on V s in cacao plantations under different levels of sun radiation. The Köhler et al. 17 modeling results were derived from a Jarvis-type sap flow model that found differences in sap flow response to changes in both vapor pressure deficit and radiation. The same occurred in López et al. 51 in Mexico, which found that variations in temperature and sun radiation significantly affected transpiration rates (and thus, sap flow rates) when there were variations in cloud cover. Studies conducted by Abdulai et al. 52 in Ghana showed the diurnal water use of the cacao crop as well as the average diurnal sap flow density, but did not specifically describe the use of models derived from environmental variables to predict sap flow. Overall, by predicting the diurnal pattern of cacao sap flow, the results from this type of modeling approach can also be used to obtain the diurnal patterns of transpiration and canopy conductance of cacao plants using the inverted Penman-Monteith equation as reported in different studies [14][15][16] . Our model determined the relation between different microclimatic conditions and cacao sap flow during the relative dry season in the Amazonian region, a period of great agronomic importance for fruit set and cacao pod production. The prediction of sap flow and water relations of cacao crops during this relatively dry season is of great economic importance. The production of cacao beans depends largely on the amount and quality of the flowers during the blooming season 53 , which in the continental Amazonian region occurs at the end of this relatively dry, hot period 35,54,55 . Cacao plants tend to be particularly active during this season and display high photosynthetic rates in order to bloom appropriately. However, unlike other cacao producing tropical regions with high diurnal insolation, the average sunshine duration in the continental Amazon is only 3-4 h per day 56 , due to its high cloud cover. This situation causes different shade configurations in cacao agroforestry systems 57 and, in turn, variations in microclimatic conditions 58 . These variations greatly affect cacao production 15 , particularly under the more shaded AF arrangements.
Our results showed that increases in PAR and VPD were positively correlated with sap flow values, while this correlation was negative between V s and RH a . This indicates that microclimatic variables decisively influenced the quantity and direction of V s during both day and night. The effect of the microclimatic conditions was modified by the upper tree canopy and plant density, which determined variations and differences in the pattern of Table 2. Mean values ± SE of the microclimatic variables during the monitoring period under different radiation intensities and at different hours of the day (official hours UTC/GMT-5:00 h): air relative humidity (RH a , %), air temperature (T a , °C), photosynthetically active radiation (PAR, µmol m −2 s −1 ), vapor pressure deficit (VPD, kPa). The agroforestry arrangements are described in the methods section and Fig. 2. Values in a column that have different letters within the same time slot indicate significant differences between AFs (posthoc LSD tests, p < 0.05).  Figure 4. Diurnal patterns of air relative humidity (RH a ); air temperature (T a ); vapor pressure deficit (VPD); and photosynthetically active radiation (PAR) under different agroforestry arrangements during the period of study. Types of AFs and legend are described in the methods section and in Fig. 2. Table 3. Pearson correlation coefficients between V s in cacao plants and microclimatic variables during the monitoring period under different radiation intensities. Types of AFs are described in the methods section and in Fig. 2. The microclimatic variables are: air relative humidity (RH a ); air temperature (T a ); photosynthetically active radiation (PAR); and vapor pressure deficit (VPD). All correlation coefficients were highly significant (p < 0.0001). www.nature.com/scientificreports/ compared to the other AFs suggests that air conditions (and probably soil humidity) were optimal at the high levels of radiation present in H PAR , which was directly related to the high demand for water of cacao trees so as to be physiologically active and perform photosynthesis 36 . In fact, the maximum photosynthetic rates in cacao trees in H PAR were significantly higher than those in cacao trees in the M PAR and L PAR AFs 57 . Besides predicting the positive sap flow in cacao plants and its differences between agroforestry arrangements, our model was also able to predict inverse or basipetal flows during the night. This inverse flow (from the main  59 and many others, who consider this phenomenon to be directly related to plant hydraulic redistribution. Hydraulic redistribution is considered indicative of a passive process of sap movement in the direction of the maximum difference in water potential between the parts of the plant (branches, trunk, roots) and the soil 59,60 . In the L PAR and M PAR arrangements, it is possible that the cacao plants had a higher water potential (closer to zero) during the night than the surface soil layers, which may have dried during the day due to water uptake by the whole plant community of the AF and soil water evaporation. Soils in these shaded AFs might have been drier than in H PAR (where there was more incident radiation, a higher temperature, and potentially more evaporation in the soil than in L PAR and M PAR ), probably due to the great demand for water Table 4. Coefficients of the best model (according with AIC and BIC criterion) predicting sap flow (V s ) as a function of a set of microclimatic variables across the different radiation intensities for cacao plantations in the Colombian Amazon. The types of AFs and legend are described in the methods section and Fig. 2. The microclimatic variables are: air relative humidity (RH a , %); air temperature (T a , °C); photosynthetically active radiation (PAR, µmol m −2 s −1 ); and vapor pressure deficit (VPD, kPa). The types of AF were included as dummy variables in the regression model to control for potentially effects of different agroforestry arrangements and incident radiation intensities on microclimatic variables and sap flow across AFs. The first three coefficients (Value) shown in the table are: the value for H PAR corresponds to β 0 (see Eq. 3) in H PAR , the value for M PAR corresponds to the difference between H PAR -β 0 and M PAR -β 0 , the value for L PAR corresponds to the difference between H PAR -β 0 and L PAR -β 0 . The following groups of coefficients (row blocks separated by a line) correspond to the coefficients for each microclimatic variable and is interpreted in the same way as described above but applied to the partial slope for each microclimatic variable. In all cases, the M PAR vs L PAR rows show the difference between their estimated coefficients.  www.nature.com/scientificreports/ by the companion plant species that make up the shade canopy in L PAR and M PAR . Our results showed that the higher the density of companion plants (L PAR > M PAR > H PAR ), the greater the hydraulic redistribution observed, suggesting there was a depletion in surface soil water content with increasing plant density. In fact, results from other studies support this idea, as soil humidity data measured in the same AFs during the dry season showed that surface soil profiles in the L PAR and M PAR AFs were, on average, drier than H PAR 61 . Moreover, the fact that our sap flow predictive models differed across AFs suggests that soil water content not only differed among AFs but also modulated cacao sap flow dynamics.

RH a T a PAR VPD
For this plant hydraulic distribution to occur, it is necessary that nocturnal transpiration stops or vastly decreases; in such cases, the highest water potential gradient in the plant soil continuum may be between the wet soil profiles and drier surface soil layers 62 and, potentially, from the trunk to the roots 63 . Low VPD values (e.g., fog and cloudy nights with a high RH a ) may impair nocturnal transpiration 62 , creating favorable conditions for water to be redistributed into dry soils. Our results align with this idea: in L PAR and M PAR , unlike in H PAR , there were low levels of VPD and very high levels of RH a , which likely inhibited the nocturnal transpiration of cacao plants, allowing the preferential sap flow to occur from the plant to the relatively drier surface soil profile. The diurnal proportion of reversed V s (i.e., the proportion of flow that went from the plant tissue to the roots and, potentially, on to the soil) was 49.3% in L PAR and 5.08% in M PAR ; these values fall within the range reported by other authors 63,64 for other forest species. This ability of cacao to redistribute water is probably relevant to reduced competition for surface water with the companion species the next day, allowing cacao trees from L PAR and M PAR to be physiologically active during the next day. Instead, and like L PAR and M PAR , the positive nocturnal sap flow exhibited by cacao trees in H PAR is probably indicative of significant plant nocturnal transpiration, which usually inhibits hydraulic redistribution 60,65 . Significant nocturnal transpiration indicates that the maximum gradient of water potential for cacao plants in H PAR would have been between the soil profile to the roots, and from the roots to the leaves and atmosphere, preventing (or reducing) hydraulic redistribution 66 .

Conclusions
The different cacao agroforestry arrangements, which differed in their structure and density of shade canopy plants, generated different incident radiation intensities, modifying the microclimatic conditions in each plantation and directly affecting diurnal patterns of sap flow in cacao trees. We accurately modelled sap flow as a function of a set of microclimatic variables, including PAR, VPD, RH a , and temperature in each AF, and accurately modelled the differences in sap flow across AFs. Nonetheless, and regardless of the type of agroforestry arrangement, the variable with the greatest effect on cacao V s was PAR (and thus, incident radiation), which also modulated RH a and T a values, having all an effect on VPD. The model can also track plant phenomena such as nocturnal transpiration (in H PAR ) and inverse nocturnal sap flow indicative of hydraulic redistribution (in M PAR and L PAR ) and, thus, could be a useful tool for managing and predicting cacao tree water use as a function of the microclimatic conditions in the different agroforestry arrangements in the Continental Amazon rainforest region.