Water stress coefficient determined by orbital remote sensing techniques

: In regions where the irrigated area is increasing and water availability is reduced, such as the West of the Bahia state, Brazil, the use of techniques that contribute to improving water use efficiency is paramount. One of the ways to improve irrigation is by improving the calculation of actual evapotranspiration (ETa), which among other factors is influenced by soil drying, so it is important to understand this relationship, which is usually accounted for in irrigation management models through the water stress coefficient (Ks). This study aimed to estimate the water stress coefficient (Ks) through information obtained via remote sensing, combined with field data. For this, a study was carried out in the municipality of São Desidério, an area located in western Bahia, using images of the Landsat-8 satellite. Ks was calculated by the relationship between crop evapotranspiration and ETa, calculated by the Simple Algorithm for Evapotranspiration Retrieving (SAFER). The Ks estimated by remote sensing showed, for the development and medium stages, average errors on the order of 5.50%. In the final stage of maize development, the errors obtained were of 23.2%.


Introduction
Agricultural activities have been responsible for the important economic growth of Brazil in recent decades (Abbade, 2014). Among the various crops that contributed to this growth, maize (Zea mays L.) is one of the main ones, being found in all regions of the country, so it is important to know the factors that influence its yield (Rodrigues et al., 2005;Rodrigues & Domingues, 2017).
Soil water availability is one of the attributes that most influence the yield of cultivated plants (Soares et al., 2011). Kresović et al. (2016) observed that water stress affected grain yield and also found the existence of a linear relationship between grain yield and actual crop evapotranspiration (ETa).
ETa is influenced, among other factors, by soil drying and can be calculated as a function of potential crop evapotranspiration (ETc) and water stress coefficient (Ks) (Rocha et al., 2014).
Several authors have studied the behavior of Ks for different crops. Rocha et al. (2014) observed that irrigation management changed according to the methodology used in the Ks calculation. Sayago et al. (2017) highlighted the feasibility of determining Ks through remote sensing images in soybean crop.
With the emergence of remote sensing and the possibility of identifying changes in surface cover through orbital sensors, potential applications in agriculture have also emerged (Cattani et al., 2017;Khanal et al., 2017;Yang et al., 2017). In irrigated agriculture, good results have been obtained in the estimation of crop coefficient by remote sensing (Alface et al., 2019;Lima et al., 2019;Sales et al., 2016Sales et al., , 2017.
Determination of reliable ETa values by means of orbital sensors makes it possible to obtain accurate Ks values. For Ks, little has been done so far in exploring the opportunities to estimate this parameter through remote sensing. Studies of this nature contribute to improving irrigation management because, in addition to providing direct information to feed the management models, they make it possible to capture the dynamics of the crop and correct the Ks values estimated by the models during crop development. Therefore, this study aimed to estimate Ks through information obtained via remote sensing, combined with field data.

Material and Methods
The study area is located in the municipality of São Desidério (Figure 1), which is in the mesoregion of the western Bahia state, belonging to the MATOPIBA (states of Maranhão, Tocantins, Piauí and Bahia, Brazil) agricultural frontier. In the center of the area, the geographic coordinates were recorded at 12º 27' 14" S and 45º 41' 16" W with altitude of 732 m. The coordinate reference system used was the Datum SIRGAS2000.
Ks was estimated in five center pivots (P1, P2, P3, P4 and P5) ( Figure 1C) cultivated with maize crop. Harvest was carried out on September 24, 2015 in all center pivots. The other pieces of information on sowing date and cycle duration are presented in Table 1. Irrigation was conducted according to the recommendation of the IRRIGER® management program (http://irriger.com.br/pt-BR).
The water stress coefficient was calculated by Eq. 1, according to Allen et al. (1998).  ETa was calculated using the SAFER (Simple Algorithm For Evapotranspiration Retrieving) model (Teixeira, 2010). For this, images of the Landsat-8 satellite, Operational Land Imager (OLI) sensor, obtained in the interval between maize sowing and harvest, were used. In this period, seven images were obtained in total, six of which were used and one discarded, due to the presence of clouds. The spatial and temporal resolutions of the images are 30 m and 16 days, respectively. The images were obtained for free on the Earth Explorer platform, by the Land Processes Distributed Active Archive (USGS, 2018).
Prior to using the images to estimate ETa, digital numbers (DN) were converted to physical values and atmospheric correction was performed concomitantly using the DOS (Dark Object Subtraction) method (Chavez, 1988). In radiometric conversion, DN is transformed into values of radiance (Eq. 2) and subsequently to reflectance at the top of the atmosphere. All processes of correction and processing of images were carried out using the computational resource QGIS 2.14.9. where: L λ -radiance at the top of the atmosphere, Wm -2 sr -1 μm -1 ; M L -band-specific multiplicative rescaling factor from the metadata file; (1) (2) A L -band-specific additive rescaling factor from the metadata (displacement); and, Q cal -spectral band pixel values (DN).
The values of gain offset are provided in the image metadata file.
The conversion from radiance to reflectance (ρλ) of the OLI instrument images (bands 1 to 7) was performed using the Eq. 3.
The Normalized Difference Vegetation Index (NDVI) was calculated by the ratio between the difference of the near infrared (NIR) and red (R) reflectances and their sum, according to Hanks (1974), using Eq. 9.
where: ρ λ -reflectance at the top of the atmosphere, dimensionless; ESUN λ -average solar irradiance at the top of the atmosphere for each band, W m -2 μm -1 ; Z -solar zenith angle (radians); and, d -Earth-Sun distance in astronomical units.
After atmospheric correction, the reflectance of the surface was obtained. The variables to estimate evapotranspiration, such as planetary albedo at the top of the atmosphere (Eq. 4) and brightness temperature (Eq. 5), were obtained adopting the procedures suggested by Teixeira et al. (2017).
where: α top -planetary albedo, dimensionless; and, ω λ -proportion of the amount of shortwave radiation from the sun at the top of the atmosphere in a particular range of the spectrum and the sum for all bands. ∑ Thermal infrared images of band 10 (spectral range from 10.6 to 11.19 μm) of the Landsat-8 satellite OLI sensor were used to estimate brightness temperature, according to Eq. 6.
Band-specific thermal conversion constants can be found in the image metadata file.
The instantaneous values of surface albedo (α 0 ) and surface temperature (T 0 ) were estimated based on the regression equations suggested by Teixeira et al. (2009) where: ρ B5 -reflectance on the wavelength intervals referring to the near infrared region, dimensionless; and, ρ B4 -reflectance on the wavelength intervals referring to the red region of the electromagnetic spectrum, dimensionless.
In the Landsat-8 satellite, these spectra refer to bands 5 and 4, respectively. The instantaneous values of the R ratio were calculated by Eq. 10.

ETa R ETo =
Potential crop evapotranspiration was calculated by Eq.12.

ETc Kc ETo =
The Kc values used to estimate ETc were obtained from Allen et al. (1998) and corrected for the conditions of the region. The values of Kc medium and Kc final were corrected using values of minimum air relative humidity (RHmin) and wind speed (U 2 ) from an automatic weather station in the study area. Data from 05/01/2009 to 04/30/2017 were used, in which the mean values of RHmin and U 2 were 46% and 1.26 m s -1 , respectively. The average height of the crop was considered to be 1.80 m in the final stage of vegetative development and 1.60 m in the final stage. Values of 1.20 and 0.60 (FAO 56) were adopted for Kc medium and Kc final , respectively.
The sensitivity of each variable of the SAFER model in Ks estimation was evaluated. For this, a multiple linear regression equation was fitted, correlating Ks with the variables NDVI, Ts and αs, according to Eq. 13. (8) where: β 0 -linear regression intercept; and, β 1 ...β n -angular coefficients linked to the variables of the SAFER model; and, X 1 ... X n -NDVI, Ts and αs variables.

Results and Discussion
The regression model, with R 2 = 0.88 and F ≤ 0.05, was significant for the influence of the variables NDVI, Ts and α S , and NDVI (P-value = 0.000) was the predictor that most influenced the Ks value, followed by Ts (P-value = 0.002) and α S (P-value = 0.004). The prediction of Ks was explained in 88% by Eq. 14.
for maize crop. Madugundu et al. (2017), in turn, found values of 0.82 in the full development stage (V12, V16 and R1) and 0.78 in the final stage (R1).
In the period of vegetative development, from 56 to 95 days after emergence (DAE), the NDVI values tended to stabilize, and image saturation may have occurred. Although very dense and high-biomass vegetation can saturate the image (Jensen, 2009), which is the case of maize crop when it reaches maximum biomass increment, Bertolin et al. (2017) stated that it was not possible to confirm this behavior when working with maize crop in the same study area, based on the scatter plot. In the present study, the same behavior was observed.
ETa ranged from 0.05 to 5.94 mm d -1 . The variation of ETa observed, especially on 06/10/2015, for P3, P4 and P5, can be attributed to the fact that maize sowing was performed on different dates. In the image of 05/25/2015, the crop was still in the process of emergence and stabilization, which justifies the absence of variation in Eta values ( Figure 2B). From 07/28/2015, it was possible to observe that the ETa values tended to become similar for the five center pivots studied. However, the last two images showed variations of 0.36 mm d -1 between P1 and P5 and 0.63 mm d -1 between P1 and the center pivots P3 and P4.
The lowest ETa values were observed at the beginning and end of the crop cycle. Low NDVI values, both at the beginning and at the end of the cycle, may explain part of this behavior, resulting from the predominance of bare soil. At the beginning of the cycle, the soil cover by the plant is small and, at the end, the leaves senescence, exposing the soil. In turn, the highest value (5.86 mm d -1 ) was found on 08/13/2015, when the crop reached maximum vegetative development. The same behavior s Ks 5.23559 1.27908NDVI 0.08950Ts 0.25696 The values of NDVI, ETa, ETc and Kc for the satellite image acquisition period are presented in Figures 2A, B, C and D, respectively.
The NDVI showed a similar response among the center pivots analyzed for most of the days studied, with values ranging from 0.24 to 0.91. The largest difference (0.21) in NDVI values was observed on 06/10, beginning of the development stage, between center pivots P1 and P5, which can be attributed to the difference in the sowing time for the five center pivots, which was equal to 7 days between P1 and P5. Similar Kc used in each center pivot was observed by Sales et al. (2016Sales et al. ( , 2017 for bean crop in the Cerrado Region of the Federal District and for tomato crop in the municipality of Silvânia, GO, Brazil, using the Landsat-8 satellite (OLI/TIRS) in both situations.
ETc ranged from 1.56 to 5.73 mm d -1 ( Figure 2C). It can be observed that the ETc value was the same in the five center pivots, on the dates of 07/27 and 08/13. This may have occurred because ETc was calculated as the product of ETo by Kc and on those dates the Kc used was the same; in the development and final stages, the Kc varied according to the DAE (0.40 to 0.51 on day 06/10; 1.05 to 1.15 on 08/29; and 0.83 to 0.93 on 09/14).
The Kc values in the stages in which it varies with time (rapid growth: 21 to 55 DAE and final: 96 to 140 DAE) were calculated using the following equations: Kc = 0.0214 DAE -0.0278, for the rapid growth stage, and Kc = -0.0136 DAE + 2.4795, for the final stage ( Figure 2D).
With the ETa values for each image and the mean daily ETc for the same date of obtaining, the water stress coefficient was calculated ( Figure 3).  The Ks values were extracted to obtain the mean values of the maps in Figure 3 and compared with the values used to perform irrigation management (Bernardo, 2019), as shown in Figure 4.
The low Ks values observed for the five center pivots in the image of the 05/25 can be explained by the low values of ETa, arising from low NDVI values in the initial period of the crop. However, at this stage the exploration by the roots occurs mainly in the surface layer. Allen et al. (1998) stated that the readily evaporable layer is up to 15 cm. Thus, it is assumed at this stage the soil could be with the dry surface, but the deeper layers, no, generating the great difference between the Ks_MD and Ks_CP.
In the image of 08/29, the Ks values estimated for the center pivots P1, P2, P3, P4 and P5 were 9, 12, 12, 12.6 and 16% lower than the Ks observed in the field, respectively. In the image of 09/14, these values were on the order of 44.6, 30.8, 33.5, 28.5 and 30.7%, respectively.
When comparing the mean ETa of the five center pivots calculated from the relationship between ETc and Ks, it was observed that the error resulting from Ks_MD originated lower ETa values for five of the six days studied, when compared with Ks_CP. On 05/25 (initial stage), the mean ETa calculated with Ks_MD was 1.46 mm d -1 lower and, on 06/10 (development), this value was 0.29 mm d -1 lower. When the plant was already in the full development stage (07/28 and 08/13), the ETa was 0.33 mm d -1 lower and 0.18 mm d -1 higher, respectively. In the final stage of development, days 08/29 and 09/14, the values were 0.57 and 1.19 mm d -1 lower, respectively.
According to Figure 4, the Ks MD showed values close to those of Ks_CP in the stages in which the crop is beginning its development (06/10) and when it is in full development (07/28 and 08/13).
This demonstrates that, although the methodology needs to be better evaluated for the initial and final stages of crop development, it has potential for use in irrigation management in the intermediate stages.
It is essential that ETa be adequately estimated, since it has a direct influence on the estimation of Ks. Any model of ETa estimation by remote sensing can be used, such as: SEBAL (Bastiaanssen et al., 1998), METRIC (Allen et al., 2007), SSEBop (Senay et al., 2013). ETa estimates by SAFER can be improved with its calibration. SAFER calibration was not performed in the present study, because the main objective was to demonstrate the potential of application of the methodology.

Conclusions
1. At the initial stage of the maize crop, the estimates of the water stress coefficient by means of remote sensing were not adequate, whereas for the development and medium stages, the estimates were close to observed value, with mean errors of 5.50%. In the final stage of maize development, the errors obtained were of 23.2%.
2. The methodology showed potential for application in irrigation management, so it needs to be better evaluated for other field situations and combined with other methods of estimating actual evapotranspiration by remote sensing.