Numerical modeling of actual evapotranspiration of a coffee crop

The evapotranspiration estimation has great importance to crop productivity and agricultural water management. In this study, evapotranspiration is analyzed in a coffee (Coffea arabica L.) crop located in Piracicaba, state of São Paulo (Brazil) using a numerical method based on the simulation of both water flow and crop activity in the unsaturated zone of the soil. Actual evapotranspiration is estimated from potential evapotranspiration using water stress functions, meteorological data, soil hydraulic parameters, crop coefficients and leaf area index values. Crop transpiration and soil evaporation are individually quantified improving the analysis of the evapotranspiration process. The numerical procedure can predict periods of crop water stress and becomes an attractive tool to analyze the effect of non-standard conditions on coffee crops and to design efficient irrigation schedules. Simulated evapotranspiration values are in good agreement with experimental values determined in the study site.


Introduction
In recent years the quantification of evapotranspiration (ET) has received special attention in many agronomic applications.ET plays a key role in the analysis of crop development under conditions of limited water resources (Rana and Katerji, 2000).The estimation of ET is also important to optimize irrigation scheduling since crop productivity depends on water availability in the soil profile during vegetative and reproductive phases.
The aim of this study is to analyze the ET in an experimental coffee crop grown in Piracicaba, state of São Paulo, Brazil, (22°43' S; 47°38' W) using a numerical method recently proposed by Cesanelli and Guarracino (2009).Several models exist in the literature to describe individual components of the land-plant-atmosphere system but unfortunately only few methods for estimating ET obtain benefit from these models.The proposed method is based on the numerical simulation of water flow in the soil profile and uses several well-known models for describing different components of this complex system.Although this method is an alternative to standard methods, it has only been tested and proven for a native green grass in La Plata (Argentina).In this context, the main contribution of this work is to extend and verify the applicability of the method to a different crop and climatic conditions, such as the case of the coffee crop grown in Piracicaba (Brazil).
Coffee (Coffea arabica L.) farming is one of the most important agribusinesses in Brazil.A number of experimental studies were conducted last years to investigate the effect of climatic and environmental conditions on coffee crop development (Carr, 2001;Villa Nova et al., 2002;DaMatta, 2004;Marin et al., 2005;Sentelhas et al., 2006;Lin, 2007;Sato et al., 2007).In this study we used numerical simulation techniques in conjunction with detailed information of soil and crop properties to evaluate the ET process.An important advantage of the proposed method is the individual computation of soil evaporation (E) and crop transpiration (T).This type of analysis is especially useful in agricultural practices where E is considered an unproductive loss of water (Wallace et al., 1999).Simulated values of ET are compared with experimental data (Silva et al., 2006) and reasonably good agreement between both series is found.The numerical results illustrate the ability of the proposed method to estimate actual ET and to predict the effective water consumption of the coffee crop.

Materials and Methods
In this section, it is presented a brief description of the numerical method and the experimental data used to estimate the E, T and ET of the coffee crop.Further details of the proposed method and its numerical implementation can be found in Cesanelli and Guarracino (2009).The experimental data are extensively described by Silva et al. (2006), Timm et al. (2006), Bruno et al. (2007), and Timm et al. (2011).
The numerical method proposed by Cesanelli and Guarracino (2009) estimates actual values of E, T and ET from their respective potential values using the following linear relationships: where subscripts A and P refer to actual and potential values, respectively; and k E , k T and k ET are coefficients that describe the effect of water stress on E, T and ET, respectively.These coefficients range from 0 to 1 and can be calculated by solving Richards' equation (Richards, 1931), which describes the water content changes in the soil profile: Sci. Agric.(Piracicaba, Braz.), v.68, n.4, p.395-399, July/August 2011 where h is the soil water matric potential head, t is the time, z is the vertical coordinate and the gravitational potential head, θ(h) is the water content, K(h) is the soil hydraulic conductivity, β T (h) is a water stress function for crop transpiration and g(z) is the root density function.Note that the righthand side of Eq. ( 2) is a sink term that describes the root water uptake by plants.
The soil evaporation process is assumed to be limited to the soil surface (z top ) and it is modeled by the following boundary condition: ( ) where P is the rainfall intensity and β E (h) is a water stress function for soil evaporation.At the bottom edge of the domain (z bot ) a Dirichlet type boundary condition is prescribed: where z wt (t) denotes the position of the water table as a function of time.
The proposed method assumes that E P , T P and ET P values are known.ET P can be expressed as (Allen et al., 1998): where: ET 0 is the reference evapotranspiration, which can be estimated from meteorological data, and k C is the crop coefficient.The potential values of E and T are obtained using the Beer-Lambert law (Huygen et al.,1997): where LAI is the leaf area index and γ is the extinction coefficient which accounts for the interception of the radiation by vegetation.
It can be shown by integrating Richards' equation in time and space that the daily water stress coefficients k E and k T have the following expressions (Cesanelli and Guarracino, 2009): where z root is the limit of the root zone and (t 0 , t 24 ) defines a time period of one day.Given that ET A is equal to the sum of E A and T A , the water stress coefficient k ET can be expressed: Based on the previous equations and definitions, the proposed method to estimate E A , T A and ET A can be stated as follows: (i) compute ET P with a daily resolution using Eq. ( 5) and obtain the partition between E P and T P using Eq. ( 6); (ii) solve Richards' equation ( 2) with the boundary conditions (3) and (4) using a finite element method; (iii) compute daily water stress coefficients k E , k T and k ET from the numerical solutions of Richards' equation using Eqs.( 7) and (8); and (iv) compute daily E A , T A and ET A using Eq. ( 1).
The experimental data used in this study were obtained in an experimental coffee plot located in Piracicaba (Silva et al. 2006).The area of the plot has 0.2-ha and is divided into 15 subplots with 120 plants each.The climate of this region is mesothermic with a dry season from April to September, being July the driest month, and a wet season from October to March.Annual average temperatures, rainfall and relative humidity are 21.1ºC, 1257 mm and 74%, respectively.The soil is a Rhodic Kandiudalf, locally called "Nitossolo Vermelho Eutroférrico" and the variety of coffee plants (Coffea arabica L.) is "Catuaí Vermelho".The period of this analysis is restricted to two years (from September 2003 to August 2005) in accordance with the available experimental data.
Reference ET (ET 0 ) is computed with the FAO Penman-Monteith equation using the standard method presented in (Allen et al., 1998) with daily meteorological data obtained from an automatic weather station located about 200 m from the experimental plot.The FAO Penman-Monteith equation is a simplified representation of the classical Penman-Monteith model that uses standard climatological records of air temperature, relative humidity, wind speed and net radiation.The range of ET 0 values varies from 0.27 to 7.35 mm per day for the study period.The daily values of ET P are calculated using Eq. ( 5) and experimental crop coefficients determined by Silva et al. (2006).The partition between E P and T P is obtained using Eq. ( 6).Due to the lack of a specific estimate of extinction coefficient γ, this parameter is assumed to be 0.53, which is the value determined by Angelocci et al. (2008) for a similar coffee crop.On the other hand, and since only few values of LAI are available, the following expression is used to estimate daily values of LAI through the year (Sumner and Jacobs, 2005): where: j is the day of year, c 1 = 2.88402, c 2 = -0.289465and c 3 = 6.6668.The values of model parameters c 1 , c 2 and c 3 are obtained by fitting expression (9) to the experimental values of LAI using an exhaustive search method.The choice of model ( 9) is motivated by the seasonal behavior of LAI data in the study region.Figure 1 shows the experimental values of LAI and the fitted curve for LAI evolution Eq. ( 9).The hydraulic properties of the soil are described using the well-known Van Genuchten model (Van Genuchten, 1980): where θ r and θ s are the residual and saturated water contents, respectively; α and n are empirical fitting parameters, Sci.Agric.(Piracicaba, Braz.), v.68, n.4, p.395-399, July/August 2011 and K S is the saturated soil hydraulic conductivity.The Van Genuchten parameters θ r , θ s , α and n are computed using the pedotransfer functions for Brazilian soils derived by Tomasella et al. (2000).Given a bulk density of 1530 kg m -3 (Timm et al., 2006) and considering the percentage of sand (29%), silt (16%), clay (55%) and organic matter (2.5%) of the Nitossolo Vermelho Eutroférrico soil (Silva et al., 2006) the following values of Van Genuchten parameters are obtained: θ s = 0.4346, θ r = 0.2275, α = 3.11 m -1 and n = 1.5492.The saturated soil hydraulic conductivity K S is estimated as 468.64 mm per day by Silva et al. (2007).
Water stress functions β E (h) and β T (h) are assumed to have the following expressions: where h E1 and h T1 are cut-off values for E and T processes; and h E2 and h T2 are cut-off values for E and T under standard conditions, respectively.Based on previous studies (e.g.Li et al., 2001) h E1 = h T1 = -150 m and h E2 = h T2 = -1 m are considered.Since more than 98% of the root system is contained in a soil layer of 1 m (Silva et al., 2006) a constant density function g(z)=1/(z top -z root ) and a root zone that extends to this depth (z root = z top -1m) are assumed.The mean water table is estimated to be 1.75 m below ground level during the simulation period.Numerical solutions of Richards' equation are obtained using a mixed finite element method in space combined with a backward Euler scheme in time.A modified Picard iteration scheme is employed to deal with the non-linear terms of Richards' equation (Celia et al., 1990).The simulation domain consists of a homogeneous soil profile of 2.0 m wide, discretized with a non-uniform mesh of 280 nodes.A variable time step is used for temporal discretization.

Results and Discussion
Daily values of both potential (T P ) and actual transpiration (T A ) for the study period are shown in Figure 2. As expected, the estimated values of T P show seasonal variations.During spring and summer T P is, on average, approximately 4 mm per day, while during autumn T P decreases and most of the values range between 1 mm per day and 3 mm per day.In particular, T P shows high variability during the summer of 2004 (values range from 1.5 mm per day to 9.5 mm per day) and an important decrease at the beginning of the winter of 2004 (values are lower than 1.5 mm per day).The temporal pattern observed for T P satisfactorily agrees with the coffee crop activity under normal conditions in Brazil, which is maximum from October to March, moderate from April to June, and minimum from July to September (Sato et al., 2007).
It is possible to identify two soil-water situations from the analysis of T A .As expected, optimal soil-water conditions are mainly observed during the wet season (October to March), when the estimates of T A and T P are similar, especially after each strong rainfall event.These optimal conditions are also observed during the first part of the dry season (April to June) because the crop transpiration demand is low (i.e.low T P ) and the water available in the soil can respond to this demand (e.g. from day 250 to day 280).Conversely, the largest differences between T A and T P take place during the last part of the dry season (July to September) indicating that water availability in the unsaturated zone is very limited (e.g. from day 370 to day 400).However, T A can abruptly increase after a strong rain event, as observed at the end of the study period (day 727).Although these specific periods of crop water stress are required by coffee plants to blossom, an extended drought can damage the plant and may cause flower loss with a significant decrease in productivity (Silva et al., 2006;Lin, 2007).The water deficit predicted by the model can be used to design irrigation management strategies so as to avoid a decline in productivity.
Simulated values of potential (E P ) and actual evaporation (E A ) are shown in Figure 3. E P values are very low compared to T P values and present a periodic pattern mainly determined by LAI evolution (Figure 1).The highest E P rates (between 0.8 mm per day and 0.9 mm per day) are observed during spring, when plants reach the minimum area covered by leaves after the period of crop water stress that takes place from July to September.From mid-summer to mid-winter, E A is approximately constant and has the lowest values (E A ≈ E P < 0.1 mm per day).During the rest of the simulation period, E A presents high variability with amplitudes lower than E P .Daily values of potential (ET P ) and actual evapotranspiration (ET A ) for the study period (Figure 4) indicate that according to the previous analysis ET A is mainly determined by crop transpiration.In spring and summer ET A presents great variability with values ranging from 1 mm per day to 8 mm per day approximately.On the other hand, during autumn and winter ET A is less variable and has very low values (ET A reaches rates higher than 3 mm per day only in a few cases).The comparison between ET A and ET P shows that water demand is satisfied over most of the simulated period, except for the last months of the dry season (July to September).Some researchers argue that, throughout these drought periods, water uptake decreases as a consequence of the small values of hydraulic conductivity in the soil profile.Note that this relation between root water uptake and soil hydraulic properties is explicitly considered in the numerical model.
Numerical estimates of ET A are validated with experimental values determined by Silva et al. (2006).The experimental data set consists of 52 ET A values obtained every 14 days from a classical water balance experiment using the available field data (measurements of rainfall, irrigation, runoff, soil water content and different meteorological variables).To compare numerical and experimental results, daily simulated values of ET A are integrated over periods of 14 days.The agreement between simulated and observed ET A values is shown in Figure 5.The root mean square difference (RMS) is 0.86 mm per day while the coefficient of determination R 2 is 0.68.Silva et al. (2006) estimated that the variability of ET values range between 27% and 470% through error propagation analysis.According to these results, simulated values can be considered within the range of the experimental determination.
Finally, it is worth mentioning that the estimation of actual ET is one of the most difficult tasks in soil and agricultural sciences due to the complex interactions amongst the components of the land-plant-atmosphere system.Although many models have been developed to characterize individual components, only a few methods for estimating ET take advantage of these models.In this sense, the proposed method integrates several of these specific models through a sink term and boundary conditions in Richards' equation.For example, the Van Genuchten model is used for describing the hydraulic properties of the soil, the FAO Penman-Monteith equation for computing reference ET, the Beer-Lambert law for partitioning E and T components, a sinusoidal function for estimating LAI evolution, and the expression (11) for modeling water stress functions.All these models are well-known and have been used in many applied studies.However, each of them can be easily replaced in the computational code by more adequate models for a particular crop, soil or climatic condition.
In conclusion, a new method for estimating ET was applied to an experimental data set of a coffee crop grown in Piracicaba.This method has been recently proposed by Cesanelli and Guarracino (2009) and is based on the numerical simulation of both water flow and crop activity in the unsaturated zone of the soil profile.Soil evaporation and crop transpiration are individually estimated improving the   understanding of the ET process.Soil hydraulic properties and crop development are included explicitly in the numerical modeling using both field data and parameter values from the literature.The proposed analysis shows that the effect of soil evaporation on ET is subtle and therefore water consumption in the experimental plot is mainly determined by crop transpiration.Two soil-water situations are identified from the comparison between potential and actual values of ET.Optimal soil-water conditions are observed during most of the simulation period while non-standard conditions are observed during specific periods of crop water stress (July to September).Simulated values of ET are in reasonable good agreement with the available experimental data.The proposed analysis illustrates the usefulness of the numerical method as a predictive and management tool in agricultural practices.