Large-Scale Spatial Modeling of Crop Coefficient and Biomass Production in Agroecosystems in Southeast Brazil

Sentinel-2 images at 10-m resolution were used for modeling crop coefficients and biomass production with the application of the so-called SAFER (Simple Algorithm for Evapotranspiration Retrieving) and Monteith model for biomass production in an area nearby the city of Águas de Santa Bárbara, in the central-western part of São Paulo State, Brazil, which presents a vast agricultural landscape mosaic, to analyze the effects of the end of the recent ENSO’s (El Niño-Southern Oscillation) most active period (2016/2017) and its posteriori effects on vegetation (until early 2018). Surface albedo, temperature, net radiation, and NDVI (Normalized Difference Vegetation Index) from the main land uses were extracted to process microclimatic comparisons. Crop coefficient (dimensionless) and biomass production (kg.ha−1.day−1) ranges for the period studied were 0.92–1.35 and 22–104 kg·ha−1·day−1 (in the area occupied by sugarcane crop), 0.56–0.94 and 15–73 kg·ha−1·day−1 (pasture), 1.17–1.56 and 25–210 kg·ha−1·day−1 (silviculture), and 1.05–1.36 and 30–134 kg·ha−1·day−1 (forest). According to the spatial and temporal consistencies, and after comparison with previous point and large-scale studies with similar climatic and thermal conditions, the SAFER and Monteith modelsshowed the ability to quantify and differentiate the large-scale crop coefficients and biomass production of different land uses in the southeast Brazil region. The SAFER algorithm with Sentinel-2 images obtained crop coefficients that indicated plant growth stages and local thermohydrological conditions at a 10-m resolution. The results are important for land use, crop yield and reforestation planning, and for water management plans for actual and future water demand scenarios, and this methodology is useful for monitoring rural and water parameters, and for precision agriculture applications.


Introduction
The multiple uses of water (agriculture, public supply, industry, etc.) and their scarcity are variables in water resource management for urban and rural areas, using the watershed as the unit for hydrological studies, as provided in Brazilian legislation.Therefore, the quantification of evapotranspiration (ET) is indispensable for water management in irrigated crops, since it represents the water demand of crops, and knowledge of it avoids both water and energy waste and productivity losses.ET is strongly influenced by vegetation type, agricultural management, environmental management, and climatic parameters [1], including solar radiation, wind, temperature, and relative humidity.
According to the 2016 Brazilian National Water Agency's (ANA) irrigation water use compilation, between 1960 and 2015, the irrigated area in Brazil increased significantly, from 462,000 hectares to 6.95 million hectares (Mha), and is predicted to expand by 45% by 2030, reaching 10 Mha [2].In Brazil, irrigated agriculture is responsible for the withdrawal of 969,000 L of water per second (969 m 3 •s −1 ), and the consumption of 745,000 L per second (745 m 3 •s −1 ); considering the other consumptive uses, these values correspond to 46% of the total withdrawal (2105 m 3 •s −1 ) and 67% of the consumption flow (1110 m 3 •s −1 ).These amounts match those of the United States where 59% of the withdrawal flow is for irrigation, and the global average of about 70% of consumption [2,3].
For these reasons, accurately estimating the amount of water withdrawals and actual consumptive use of water for irrigation at a regional scale is necessary; however, it is a difficult task.Usually, the amount of irrigation water is calculated based on actual evapotranspiration (ET a ) using a reference evapotranspiration (ET 0 ) multiplied by a crop coefficient (Kc) [1].Kc curves have been established according to their correspondence to the growth stage of the crop [1,4], leaf area index (LAI) [5], accumulated days-degrees (DDac) [6,7], and as a function of time [8].However, this method, because of its simplicity, disregards the influence of spatial factors in determining ET a .Therefore, obtaining ET by satellite imagery providing spatially distributed information at low cost is promising [9].
Over the last two decades, remote sensing techniques have enabled monitoring at different spatial and temporal scales in the various biomes of the world to determine the biophysical parameters at the surface and the ET a at different spatial and temporal scales, for example, SEBAL (Surface Energy Balance Algorithm for Land), elaborated by Bastiaanssen et al. [10]; S-SEBI (Surface Energy Balance Index) by Roerink et al. [11]; and SEBS (Surface Energy Balance System) by Su [12].A combination of these methodologies is SAFER (Simple Algorithm for Evapotranspiration Retrieving) developed by Teixeira et al. [13], which models the ratio between actual (ET a ) and reference (ET 0 ) evapotranspiration, which is equivalent to the crop coefficient (Kc), or the so-called evaporative fraction (ET r ) by some authors [6,7,9,13,14].
SAFER is a simplified algorithm that has shown good results in the estimation of large-scale ET.Hernandez et al. [14] compared the SAFER and SEBAL algorithms to apply the crop coefficient approach [1] in the northwestern part of São Paulo State, Brazil on corn, beans, and sugarcane crops irrigated by center pivots.In that study, SAFER presented better performance than SEBAL under the conditions of low soil cover by the canopies and similar performance under high soil cover by canopies, which indicates that the SAFER model is better suited to monitor crops from early (with lower soil cover) to near-harvest stages (higher soil cover), as has been reported [6,7].This algorithm has the advantage of avoiding the use of the thermal band, requiring only the measured reference evapotranspiration (ET 0 ), global radiation (R G ), and average air temperature (T a ), and radiometric data from remote sensing tools, to reach the energy balance including the ET a , which can be used as a tool to help irrigation management in different crops by spatially defining the ET r ratio or crop coefficient.
Linked to ET, there is the production of biomass, which is fundamental for understanding the dynamics of energy and carbon fluxes in soil [15] in agroecosystems with varied land uses.The relationship between radiometric data obtained from a canopy or vegetated field and parameters that categorize the growth stages of plants is the principle used for the study of biomass production by remote sensing.The combination of measurements of electromagnetic radiation reflected by vegetation in some spectral bands results in vegetation indices (VI).Using vegetation indexes, it is possible to reduce the dimensionality of the data and to highlight the spectral response of the vegetation when related to the soil [16].Additionally, it is necessary to study the production of biomass in an agroecosystem context, where several uses of the land interact.For Rezende et al. [17], the scarcity of information on the primary production of biomass in the Cerrado is mainly related to the great diversity of species and also to the high variability of the shape of the tree trunk and crown, even of trees of the same species.The model of biomass production proposed by Monteith [18], from the quantification of solar radiation, can be applied in conjunction with satellite data for predictions of biomass growth and crop productivity [18][19][20][21].
The availability of freely moderate-to-high spatial resolution satellite imagery stimulated the use of remote sensing solutions in spatial modeling, especially after the launch of the Landsat-8 satellite in 2013 [22].While Landsat-8 imagery (at 30 m spatial resolution) has been widely used for ET mapping [6,7, [23][24][25], other remote sensing imagery, like RapidEye [26] and MODIS [27], with higher and lower spatial resolution than the Landsat-8, respectively, began to be used successfully.A new option is the Sentinel-2 satellite, launched in June 2015, which provides multi-spectral images in 13 bands at different spatial resolutions, including three visible (RGB) and one near infrared (NIR) bands at 10 m spatial resolution, six red-edge and shortwave infrared (SWIR) bands at 20 m, and three atmospheric correction bands at 60 m [28] acquired at a swath width of 290 km and a five-day temporal resolution.
The Sentinel-2 and Landsat-8 satellites have similar band configurations, highlighting that the bands referring to the visible and the near infrared have near wavelength ends in the two satellites (Figure 1), which allows the use of both for the same objectives [22,28].However, Sentinel-2 gathers a higher spatial and temporal resolution, which is an advantage.
Horticulturae 2018, 4, x 3 of 20 even of trees of the same species.The model of biomass production proposed by Monteith [18], from the quantification of solar radiation, can be applied in conjunction with satellite data for predictions of biomass growth and crop productivity [18][19][20][21].
The availability of freely moderate-to-high spatial resolution satellite imagery stimulated the use of remote sensing solutions in spatial modeling, especially after the launch of the Landsat-8 satellite in 2013 [22].While Landsat-8 imagery (at 30 m spatial resolution) has been widely used for ET mapping [6,7, [23][24][25], other remote sensing imagery, like RapidEye [26] and MODIS [27], with higher and lower spatial resolution than the Landsat-8, respectively, began to be used successfully.A new option is the Sentinel-2 satellite, launched in June 2015, which provides multi-spectral images in 13 bands at different spatial resolutions, including three visible (RGB) and one near infrared (NIR) bands at 10 m spatial resolution, six red-edge and shortwave infrared (SWIR) bands at 20 m, and three atmospheric correction bands at 60 m [28] acquired at a swath width of 290 km and a five-day temporal resolution.
The Sentinel-2 and Landsat-8 satellites have similar band configurations, highlighting that the bands referring to the visible and the near infrared have near wavelength ends in the two satellites (Figure 1), which allows the use of both for the same objectives [22,28].However, Sentinel-2 gathers a higher spatial and temporal resolution, which is an advantage.

Data Used
The dates of cloudless images (and their corresponding day of year, DOY) were: 16 (70), and 16 March (75).In 2018, the temporal resolution of the Sentinel-2 satellite went from 15 to 5 days.Bands in the visible region of the electromagnetic spectrum (2-4) and infrared (8) were used.Digital image processing was done using the software, R 3.5.0[29], a statistical software composed of a language and an environment focused on numerical and statistical The main goals of the current research were to apply the SAFER and Monteith modelsusing Sentinel-2 images together with agrometeorological data, covering different seasons from March 2017 to March 2018.Thus, this study aimed to quantify the large-scale crop coefficient and biomass production in an area near the city of Águas de Santa Bárbara, in the central-western part of São Paulo State, Brazil (which presents a vast agricultural landscape mosaic) during the end of the recent ENSO's (El Niño-Southern Oscillation) most active period (2016/2017) and verify its posteriori effects on vegetation (until early 2018).

Data Used
The dates of cloudless images (and their corresponding day of year, DOY) were: 16 (70), and 16 March (75).In 2018, the temporal resolution of the Sentinel-2 satellite went from 15 to 5 days.Bands in the visible region of the electromagnetic spectrum (2-4) and infrared (8) were used.Digital image processing was done using the software, R 3.5.0[29], a statistical software composed of a language and an environment focused on numerical and statistical calculations.The packages used were 'raster', 'sp', and 'rgdal'.Final maps were made using the GIS-software, ArcMap™ 10.4.1 [30].
For each DOY, land use classification was performed to make zonal statistics consistent with land use dynamics.The technique used was the Multivariate Maximum Likelihood Classification to identify the thematic classes.This technique uses the mean and covariance of samples to calculate the statistical probability of an unknown pixel belonging to a given field-recognized class.After this probabilistic evaluation, the pixel is assigned to the category with the highest probability of compatibility [31].

Study Area
Watersheds located at 24 • 48 S, 49 • 13 E, in the municipality of Águas de Santa Bárbara, in the central-western part of São Paulo State, Brazil, were studied.The agrometeorological station within this area has a reference surface where the albedo is 0.23 and the LAI is 2.88.The land use of these watersheds in this research (based on DOY 75/2017), which is presented in Figure 2, are varied: natural vegetation (Cerrado biome) (67 km 2 ), silviculture (Eucalyptus) (22 km 2 ), sugarcane crops (25 km 2 ), urban area (2 km 2 ), regenerating fields (5 km 2 ), roads (0.4 km 2 ), and pasture (40 km 2 ).The northern region of the study area is largely occupied by silviculture.The western region is occupied by agriculture (sugarcane).During the year, there is harvesting of Eucalyptus and sugarcane in these areas, increasing the area of bare soil.In the southern region, there is an urban area and the pasture and forest areas are in the central and eastern region of the study area.Considering this dynamic spatial configuration of land uses, remote sensing tools are suitable for evaluating the water and vegetation conditions for developing policies for agricultural water management.
For each DOY, land use classification was performed to make zonal statistics consistent with land use dynamics.The technique used was the Multivariate Maximum Likelihood Classification to identify the thematic classes.This technique uses the mean and covariance of samples to calculate the statistical probability of an unknown pixel belonging to a given field-recognized class.After this probabilistic evaluation, the pixel is assigned to the category with the highest probability of compatibility [31].

Study Area
Watersheds located at 24°48′ S, 49°13′ E, in the municipality of Águas de Santa Bárbara, in the central-western part of São Paulo State, Brazil, were studied.The agrometeorological station within this area has a reference surface where the albedo is 0.23 and the LAI is 2.88.The land use of these watersheds in this research (based on DOY 75/2017), which is presented in Figure 2, are varied: natural vegetation (Cerrado biome) (67 km 2 ), silviculture (Eucalyptus) (22 km 2 ), sugarcane crops (25 km 2 ), urban area (2 km 2 ), regenerating fields (5 km 2 ), roads (0.4 km 2 ), and pasture (40 km 2 ).The northern region of the study area is largely occupied by silviculture.The western region is occupied by agriculture (sugarcane).During the year, there is harvesting of Eucalyptus and sugarcane in these areas, increasing the area of bare soil.In the southern region, there is an urban area and the pasture and forest areas are in the central and eastern region of the study area.Considering this dynamic spatial configuration of land uses, remote sensing tools are suitable for evaluating the water and vegetation conditions for developing policies for agricultural water management.According to the Köppen classification, the climate of the Águas de Santa Bárbara region is tropical subhumid (Cwa-warm climate with dry winter).For all months of the year (except July), the monthly average temperature exceeds 18 °C, with the warmest and coldest months being February (average of 24.6° C) and July (average of 17.9 °C), respectively.The average total annual rainfall is 1353 mm, with the wettest and driest months being January (average of 196.4 mm) and July (average of 41.4 mm) [32], respectively.
Spatially understanding large-scale microclimatic relationships between silviculture (Eucalyptus), pasture, natural vegetation, and sugarcane crops is important for the rational use of natural resources, and the consequent reduction of negative environmental impacts, supporting policy planning and decision-making on natural resources.The difficulties of measuring and analyzing these components only by field measurements have highlighted the importance of remote sensing coupling and agrometeorological data.According to the Köppen classification, the climate of the Águas de Santa Bárbara region is tropical subhumid (Cwa-warm climate with dry winter).For all months of the year (except July), the monthly average temperature exceeds 18 • C, with the warmest and coldest months being February (average of 24.6 • C) and July (average of 17.9 • C), respectively.The average total annual rainfall is 1353 mm, with the wettest and driest months being January (average of 196.4 mm) and July (average of 41.4 mm) [32], respectively.
Spatially understanding large-scale microclimatic relationships between silviculture (Eucalyptus), pasture, natural vegetation, and sugarcane crops is important for the rational use of natural resources, and the consequent reduction of negative environmental impacts, supporting policy planning and decision-making on natural resources.The difficulties of measuring and analyzing these components only by field measurements have highlighted the importance of remote sensing coupling and agrometeorological data.

Simple Algorithm for Evapotranspiration Retrivieng (SAFER)
All the regression coefficients of the parameters in Figure 3 and equations below were determined in different areas of Brazil with digital images from different satellites, like Landsat [6,7,13,14,23-25], RapidEye [26], and MODIS [9,27], with field measurements for calibration involving strong contrasting agroecosystems, and under different thermohydrological conditions throughout several years as semiarid areas [9,13,25], wetlands [27], tropical [23], and subtropical areas [6,7,14].Most of the regressions performed from Landsat-8 were used in this work due to their similarity to the radiometric data on the visible (RGB) and near infrared (NIR) range of the Sentinel-2 satellite [22,28].

Simple Algorithm for Evapotranspiration Retrivieng (SAFER)
All the regression coefficients of the parameters in Figure 3 and equations below were determined in different areas of Brazil with digital images from different satellites, like Landsat [6,7,13,14,23-25], RapidEye [26], and MODIS [9,27], with field measurements for calibration involving strong contrasting agroecosystems, and under different thermohydrological conditions throughout several years as semiarid areas [9,13,25], wetlands [27], tropical [23], and subtropical areas [6,7,14].Most of the regressions performed from Landsat-8 were used in this work due to their similarity to the radiometric data on the visible (RGB) and near infrared (NIR) range of the Sentinel-2 satellite [22,28].Albedo is defined as the ratio between reflected and incident sunlight, and is an important parameter in the study of climate change, desertification, fires, and environmental impacts [33].For retrieving surface albedo ( ), as shown in Figure 2, firstly, the planetary albedo for the entire solar spectrum ( ) was calculated as the total sum of the different narrow-band reflectances ( ) values according to weights for each band ( ), according to Equation (1): The weights for the different bands were computed as the ration of the amount of the incoming shortwave radiation from the sum in a particular band and the sum of incoming shortwave radiation for all the bands at the top of the atmosphere (TOA).The daily  values were obtained according to Equation (2): where  and  are regression coefficients, which, for a 24-h period, were considered as 1.70 and 0.13, obtained from field and satellite measurements [6, 34,35].
The NDVI (Normalized Difference Vegetation Index) is an indicator related to the land cover [36] obtained from a satellite image according to Equation (3): Albedo is defined as the ratio between reflected and incident sunlight, and is an important parameter in the study of climate change, desertification, fires, and environmental impacts [33].For retrieving surface albedo (α 0 ), as shown in Figure 2, firstly, the planetary albedo for the entire solar spectrum (α P ) was calculated as the total sum of the different narrow-band reflectances (r band ) values according to weights for each band (w band ), according to Equation (1): The weights for the different bands were computed as the ration of the amount of the incoming shortwave radiation from the sum in a particular band and the sum of incoming shortwave radiation for all the bands at the top of the atmosphere (TOA).The daily α 0 values were obtained according to Equation (2): where a and b are regression coefficients, which, for a 24-h period, were considered as 1.70 and 0.13, obtained from field and satellite measurements [6, 34,35].
The NDVI (Normalized Difference Vegetation Index) is an indicator related to the land cover [36] obtained from a satellite image according to Equation (3): where r N IR and r RED represent the reflectance over the range of wavelengths in the near infrared (NIR) and red (RED) regions of the solar spectrum, respectively.Net radiation (R N , W•m −2.sr −1.µm −1 ) was obtained by Slob's equation [23,37] in Equation ( 4), which was derived by Teixeira et al. [34], considering four energy balance field experiments involving different thermohydrological conditions in Brazil: where τ sw is the shortwave atmosphere transmissivity defined as 44% of R G to the incident solar radiation at the top of atmosphere [27] and the a L coefficient from Equation ( 5) can be explained by variations in 24-h T A values [6] by Equation ( 8): where T A is the average air temperature, R G is the 24-h values of global solar radiation, and c and d are regression coefficients equal to 6.99 and 39.93 [6], respectively.Bastiaanssen et al. [10] applied a constant value of a L = 110 without considering the local thermal conditions.The residual method for retrieving the surface temperature (T S ) was proposed by Teixeira et al. [27] according to Equation (6): where R G is the solar radiation, R N is the net radiation, T A is the average air temperature, ε A and ε S are the atmospheric and surface emissivities, respectively, and σ is the Stefan-Boltzmann constant (5.67 × 10 −8 Wm −2 K −4 ).The denominator of the fraction of Equation ( 6) is a model of the long-wave atmospheric radiation [23].We reduced the equation by applying Equation (4) in Equation (6), resulting in Equation ( 7), which depends only on the second term of Equation ( 4), without an uncertainty increase: where ε A and ε S are the atmospheric and surface emissivities, respectively, and σ is the Stefan-Boltzmann constant (5.67 × 10 −8 Wm −2 K −4 ).Emissivity is the ability of a surface to emit absorbed energy [13].A dark land surface absorbs more R G , and has a higher R N than a bright one.Teixeira et al. [6,37] performed regression analyzes for atmospheric (ε A ) and surface (ε S ) emissivity through transmissivity and NDVI values in Equations ( 8) and ( 9), respectively: where a A , b A , a S , and b S are regressions, which, from Teixeira et al. [34], were considered, respectively, as 0.94, 0.10, 0.06, and 1.00.The application of this methodology for surface temperature retrieval provides 10-m resolution results, while the SEBAL requires the use of thermal bands, which in Sentinel-2, has a spatial resolution of 1 km.The net radiation was also evaluated in this work as an indicator of the energy supply to the soil [25].It was decided that the residual method would be used to obtain the surface temperature to preserve the spatial resolution, since the thermal bands have a resolution of 1 km, which increases the uncertainty in the modeling [38].
The crop coefficient (Kc) or the ratio between actual ET (ET a ) and potential ET (ET 0 ), the so-called evaporative fraction (ET r ), was calculated according to Equation (10): where e and f are regression coefficients, 1.8 and −0.008, respectively.An advantage of Equation ( 9), compared with Bastiaanssen et al. [10], is the possibility of application in any environment according to the thermal conditions without the need for an anchor pixel (so-called "hot pixel" and "cold pixel" in the SEBAL model, which requires two points with extreme and opposite thermal conditions in the same analysis, so SAFER, different from SEBAL, allows good analysis in totally irrigated or totally dry locations [7,9,13]).Teixeira [39] has demonstrated that there are no significant differences between the daily and satellite overpass values of the ET r ratio.The average ET r values in sugarcane irrigated areas gave Kc, allowing modelling of their values with DDac (accumulated days-degrees, • C•day −1 ) using the basal temperature of 18 • C [40] and a regression polynomial equation determined by the joint use of agrometeorological and satellite data [6].According to Teixeira [13], the advantage of the use of this relationship is the transfer of Kc values to different thermal conditions.

Monteith Model of Biomass Production
The Absorbed Photosynthetically Active Radiation (APAR) was approximated directly from PAR according to Equation (11): where f PAR was estimated from the NDVI values [13,41] according to Equation (12): where the coefficients, i and j, were considered, respectively, to be 1.257 and −0.161 [6, 14,24,26,41], and the BIO was quantified according to Equation (13) [6, 14,24]: where ε max is the maximum light use efficiency, which was considered as 2.5 g•MJ −1 for the majority of C4 species [35], and 0.864 is a unit conversion factor [13].

Weather Drivers
Weather conditions influence the photosynthetic activity of crops and water fluxes.For this reason, soil moisture is essential to maintain yields at optimum levels, as water stress can harm vegetative growth and crop development of Eucalyptus (silviculture) [42], sugarcane [43], and pasture [44].Solar radiation (R G ), considered the directed and undirected radiation from the sun, was integrated over all wavelengths in the shortwave interval.Rainfall considered the total daily precipitation.Average temperature considered the mean value of air temperature at a 24 h scale.
Figure 4 shows the rainfall, average temperature, and solar radiation in the studied period.The variable with the greatest variation was precipitation.Between March and June 2017, there were high-intensity precipitations (with rainfall over 100 mm in March and over 50 mm between April and June) followed by long dry periods and solar radiation levels above 30 MJ•day −1 , while between July 2017 and April 2018, there were lower intensity precipitations and less rainfall spacing with lower levels of solar radiation (between 8 and 20 MJ•day −1 ).The delay and concentration of rains in 2017 are effects of the El Niño-Southern Oscillation [45].DOY 75, 135, and 175/2017 occurred in a period of intense precipitation and high solar radiation; DOY 235 and 270/2017 followed days of intense rainfall and low solar radiation, while the DOY of 2018 followed days of low intensity rain, more uniform distribution, and less incident solar radiation.

Input Parameters of SAFER Algorithm
The spatial results are presented in the form of maps, and in mean values with standard deviation for each land use.In this way, it is possible to analyze spatially and numerically, thus making the analysis more easily assimilated and self-connected.The areas of sugarcane, pasture, silviculture, and forest were considered as relevant land uses.Areas of exposed soil (mainly areas of sugar cane and silviculture after harvest), fields in regeneration, roads, and water bodies were not accounted for.The results have a spatial resolution of 10 m.
Figure 5 shows the spatial distribution of the surface albedo for each DOY and Table 1 presents the mean values and standard deviation for each DOY and relevant land use.

Input Parameters of SAFER Algorithm
The spatial results are presented in the form of maps, and in mean values with standard deviation for each land use.In this way, it is possible to analyze spatially and numerically, thus making the analysis more easily assimilated and self-connected.The areas of sugarcane, pasture, silviculture, and forest were considered as relevant land uses.Areas of exposed soil (mainly areas of sugar cane and silviculture after harvest), fields in regeneration, roads, and water bodies were not accounted for.The results have a spatial resolution of 10 m.
Figure 5 shows the spatial distribution of the surface albedo for each DOY and Table 1 presents the mean values and standard deviation for each DOY and relevant land use.making the analysis more easily assimilated and self-connected.The areas of sugarcane, pasture, silviculture, and forest were considered as relevant land uses.Areas of exposed soil (mainly areas of sugar cane and silviculture after harvest), fields in regeneration, roads, and water bodies were not accounted for.The results have a spatial resolution of 10 m.
Figure 5 shows the spatial distribution of the surface albedo for each DOY and Table 1 presents the mean values and standard deviation for each DOY and relevant land use.The albedo presented low spatial variability and was associated with land use and the heterogeneous vegetation cover, while the temporal variation was affected by the local climatic changes besides the presence of the anthropic action by the management of the adjacent agricultural areas.The albedo was larger and had less temporal variation in the pasture, while the lowest values were found in the forest, also with a low temporal variation.In the areas of sugar cane and silviculture, the tendency of albedo values was to increase with the development of the crop due to the increase of the soil cover, being evident in the relation between the albedo and the leaf area as reported by André et al. [46], which found the average albedo value for the whole cycle of the sugarcane crop as 0.20 ± 0.029.Giongo and Vettorazzi [47] diagnosed higher albedo in pasture areas, compared to sugarcane, due to its smaller leaf area.Giongo et al. [48] obtained the albedo in a sugarcane area, with values ranging from 0.13 to 0.23.Cabral et al. [49] reported values between 0.12 and 0.32 for non-irrigated sugarcane depending on the growth stage of the crop, corroborating results obtained in this research for areas of sugarcane.The area with Eucalyptus presented albedo values between 0.16 and 0.19, a low amplitude of the values.This reflects the small modification in characteristics because it is a perennial system of culture and can be compared to areas of forest or dense vegetation, which have presented similar behavior, registering values between 0.12 and 0.17.Boegh et al. [50], in Denmark, using TM-Landsat 5 images in areas of dense vegetation in summer, when the climate in Denmark is similar with autumn and spring in Brazil, a value of 0.18 for albedo measurement was obtained.Another cause for seasonal variations of surface albedo is the surface moisture during the rainy season, when the leaves were alternatively wet and dry due to water interception.Many authors [6, 34,35,[51][52][53] have found linear relationships between albedo and surface moisture.
Figure 6 shows the spatial distribution of the net radiation for each DOY and Table 2 presents the mean values and standard deviation for each DOY and relevant land use.Net radiation is strongly dependent on solar radiation (RG), and the capacity to differentiate between land uses through this variable is low, as already diagnosed by Teixeira et al. [25,27].However, there was a tendency for a higher radiation balance in the areas occupied by Eucalyptus and forest and lower values in pasture and sugarcane.Menezes et al. [54] found the same behavior and obtained values of radiation balance for areas of Eucalyptus and forest in Southeastern Brazil that corroborate this tendency and also the dependence of net radiation with incident solar radiation, since the authors measured a value corresponding to twice the solar radiation of our study and also reported double the balance of radiation obtained here.Gomes et al. [55] found average values of net radiation for sugarcane crops, native vegetation (Cerrado), and Eucalyptus of 10.3, 11.4, and 10.8 MJ.day −1 , respectively.
Figure 7 shows the spatial distribution of the surface temperature for each DOY and Table 3 presents the mean values and standard deviation for each DOY and relevant land use.
The use of the residual method for surface temperature retrieval (Equation ( 7)), without the use of the thermal bands currently available by the Sentinel-3 satellite that has a spatial resolution of 1 km, allowed the spatial analysis to be performed with 10 m of spatial resolution as well as the other input parameters.Silviculture and forest presented surface temperatures around 1 K smaller than pasture and sugarcane, agreeing with Meneses et al. [54], which also verified the decrease of surface temperature with the expansion of Eucalyptus.Surface temperature is mainly dependent on the absorbed solar radiation, which is converted into heat energy by the transfer of long-wave radiation from their surfaces to the lower atmosphere.Net radiation is strongly dependent on solar radiation (R G ), and the capacity to differentiate between land uses through this variable is low, as already diagnosed by Teixeira et al. [25,27].However, there was a tendency for a higher radiation balance in the areas occupied by Eucalyptus and forest and lower values in pasture and sugarcane.Menezes et al. [54] found the same behavior and obtained values of radiation balance for areas of Eucalyptus and forest in Southeastern Brazil that corroborate this tendency and also the dependence of net radiation with incident solar radiation, since the authors measured a value corresponding to twice the solar radiation of our study and also reported double the balance of radiation obtained here.Gomes et al. [55] found average values of net radiation for sugarcane crops, native vegetation (Cerrado), and Eucalyptus of 10.3, 11.4, and 10.8 MJ•day −1 , respectively.
Figure 7 shows the spatial distribution of the surface temperature for each DOY and Table 3 presents the mean values and standard deviation for each DOY and relevant land use.  Figure 8 shows the spatial distribution of the NDVI for each DOY and Table 4 presents the mean values and standard deviation for each DOY and relevant land use.The use of the residual method for surface temperature retrieval (Equation ( 7)), without the use of the thermal bands currently available by the Sentinel-3 satellite that has a spatial resolution of 1 km, allowed the spatial analysis to be performed with 10 m of spatial resolution as well as the other input parameters.Silviculture and forest presented surface temperatures around 1 K smaller than pasture and sugarcane, agreeing with Meneses et al. [54], which also verified the decrease of surface temperature with the expansion of Eucalyptus.Surface temperature is mainly dependent on the absorbed solar radiation, which is converted into heat energy by the transfer of long-wave radiation from their surfaces to the lower atmosphere.
Figure 8 shows the spatial distribution of the NDVI for each DOY and Table 4 presents the mean values and standard deviation for each DOY and relevant land use.Other authors obtained similar NDVI values in southeast Brazil, such as Fernandes et al. [43], who found average values of NDVI for sugarcane between 0.4 and 0.7; Pereira et al. [56], who found average values of NDVI for sugarcane between 0.4 and 0.75; and Castanheira et al. [57], who found an average value of 0.67 for silviculture.Data normalization has made vegetation features and characteristics more discernible in comparison to other input parameters.After rainy days with low levels of solar radiation, NDVI values decreased in all land uses compared to dry days with high levels of solar radiation, as was also reported in Lucas and Schuler [58] and Pereira et al. [56].Thus, by analyzing the surface albedo, NDVI, and surface temperature simultaneously, spatial variations  Other authors obtained similar NDVI values in southeast Brazil, such as Fernandes et al. [43], who found average values of NDVI for sugarcane between 0.4 and 0.7; Pereira et al. [56], who found average values of NDVI for sugarcane between 0.4 and 0.75; and Castanheira et al. [57], who found an average value of 0.67 for silviculture.Data normalization has made vegetation features and characteristics more discernible in comparison to other input parameters.After rainy days with low levels of solar radiation, NDVI values decreased in all land uses compared to dry days with high levels of solar radiation, as was also reported in Lucas and Schuler [58] and Pereira et al. [56].Thus, by analyzing the surface albedo, NDVI, and surface temperature simultaneously, spatial variations can be mainly attributed to variations in R G and soil moisture [6].

Crop Coefficient
Figure 9 shows the spatial distribution of the crop coefficient for each DOY and Table 5 presents the mean values and standard deviation for each DOY and relevant land use.

Crop Coefficient
Figure 9 shows the spatial distribution of the crop coefficient for each DOY and Table 5 presents the mean values and standard deviation for each DOY and relevant land use.DOYs 135, 175, and 235 presented the highest ET r values because they are after intense precipitations and occur during the phenological stages of sugarcane and Eucalyptus with the highest water consumption.ET r characterizes the water status in the root zones [59].The most important reason for the highest ET r values were prior rainy days, which make the soil more moist, but it also depends on stomatal regulation and plant adaptation to water scarcity conditions [60].Zhou and Zhou [61] concluded that air humidity, air temperature, and the available energy were the most important variables for the Kc variations in a reed marsh located in Northeast China.The region in the north of the study area occupied by silviculture reaches the highest ET r values in DOY 175, before the beginning of the harvest season.The sugarcane areas located in the west of the study area have their highest ET r in DOY 175, before harvest, when ET r decreases in this region, as seen in DOY 235.Forest and pasture areas show ET r stability.
The crop coefficient of Eucalyptus varied between 1.17 and 1.56, following the crop coefficient presented by Queiroz et al. [42].The crop coefficient of the pasture varied between 0.56 and 0.81, close to that of Capim-Elefante (Pennisetum purpureum) as described by Muniz et al. [62], which obtained an average Kc between 0.45 and 0.78, and also close to Capim-Tanzânia (Panicum maximum) as described by Bueno et al. (2009), which obtained an average Kc between 0.5 and 0.98 in the southeast Brazil region.Lima et al. [63] found crop coefficient values for forest ("Cerrado") between 1.00 and 1.10.
Figure 10 shows the relationship between mean values of crop coefficient (ET r ) and accumulated degree-days (DDac) in the sugarcane crops in the studied period.The crop coefficient of Eucalyptus varied between 1.17 and 1.56, following the crop coefficient presented by Queiroz et al. [42].The crop coefficient of the pasture varied between 0.56 and 0.81, close to that of Capim-Elefante (Pennisetum purpureum) as described by Muniz et al. [62], which obtained an average Kc between 0.45 and 0.78, and also close to Capim-Tanzânia (Panicum maximum) as described by Bueno et al. (2009), which obtained an average Kc between 0.5 and 0.98 in the southeast Brazil region.Lima et al. [63] found crop coefficient values for forest ("Cerrado") between 1.00 and 1.10.
Figure 10 shows the relationship between mean values of crop coefficient (ETr) and accumulated degree-days (DDac) in the sugarcane crops in the studied period.ETr is a measurement of the soil moisture conditions and can indicate the impact of a water deficit on the crop.Its trend was also quantified throughout a generalized growing season as a function of DDac (Figure 8).ETr at different sugarcane crop stages was between 0.2 and 1.25.This range agrees with the values in the standard work of FAO (Food and Agriculture Organization) [1], which justified the confidence in using the SAFER algorithm for the Kc approach.Also, it could be concluded that in most of the critical periods of water use during the sugarcane crop stages, between DDac of 750 and 1750 °C•day −1 , the ETr values were above 0.80 and even higher than 1.00 because of the additional energy supply by horizontal heat advection, as indicated by Teixeira et al. [6,35].ET r is a measurement of the soil moisture conditions and can indicate the impact of a water deficit on the crop.Its trend was also quantified throughout a generalized growing season as a function of DDac (Figure 8).ET r at different sugarcane crop stages was between 0.2 and 1.25.This range agrees with the values in the standard work of FAO (Food and Agriculture Organization) [1], which justified the confidence in using the SAFER algorithm for the Kc approach.Also, it could be concluded that in most of the critical periods of water use during the sugarcane crop stages, between DDac of 750 and 1750 • C•day −1 , the ET r values were above 0.80 and even higher than 1.00 because of the additional energy supply by horizontal heat advection, as indicated by Teixeira et al. [6,35].

Biomass Production
Figure 11 shows the spatial distribution of the biomass production (BIO) for each DOY and Table 6 presents the mean values and standard deviation for each DOY and relevant land use.

Biomass Production
Figure 11 shows the spatial distribution of the biomass production (BIO) for each DOY and Table 6 presents the mean values and standard deviation for each DOY and relevant land use.The standard deviation was high, indicating low uniformity.The biomass production was dependent on RG, local vegetation, and water available in the soil [27], and there was a relationship  The standard deviation was high, indicating low uniformity.The biomass production was dependent on R G , local vegetation, and water available in the soil [27], and there was a relationship between ET r and BIO in Equation ( 13); the periods with the highest BIO were the same as those for ET r .Santana et al. [64] found evidence of the correlation of production in coniferous forests with available water in the soil.High rainfall providing good moisture levels in the root zone, together with large R G levels, increase photosynthetic activity and favors large biomass production rates.Between days 235 and 270, there was a decrease in ET r and NDVI values due to the withdrawal and, consequently, a reduction of areas occupied by silviculture and sugarcane, this being evident in DOY 270, in which there were no pixels with a value greater than 45.2 kg•ha −1 •day −1 .
Oliver and Singels [65] reported strong BIO declines on sugarcane crops under irrigation conditions, when water application was reduced by around 50%.Andrade et al. [66] obtained average BIO of sugarcane crops values between 100 and 160 kg•ha −1 •day −1 during several crop periods and also affirmed that the soil moisture effects on BIO vary according to planting and harvesting dates.In South Africa, Donaldson et al. [67] reported seasonal variations affecting BIO in a number of sugarcane cultivars in the ranges between 90 and 184 kg•ha −1 •day −1 .
BIO is directly related to the product harvested through the harvest index (HI), the percentage of biomass effectively harvested and transformed into the final product [6,14,24-27], and the fraction of carbon in the biomass [25][26][27].The HI depends on whether the property is mechanized or not, and economic and topographical factors [9].ET r and BIO modeling, besides comparing different land uses as made in this paper, can compare different pasture managements [68], the effect of irrigation on crops [9,14,23], and the impact of drought on large areas [69].ET r and BIO at a 10-m resolution are suitable for precision agriculture applications [70,71].
Despite not having simultaneous field measurements for validations during the years of 2017 and 2018, because the study area was not yet equipped with lysimeters or eddy covariance, the similarities of the current results with those from the literature provide confidence for the application of the joint use of the SAFER algorithm and Monteith model to Sentinel-2 images without the thermal band in southeast Brazil, assuring a spatial resolution of 10 m for all components, showing that the calibrations made in other papers present consistent results for the study area.In future work, new calibrations and validations can be performed using simultaneous eddy covariance data for actual evapotranspiration (following the regression methodology of Teixeira [39]) or lysimeters to calibrate Equation ( 9) to the study area, and also other regression equations when needed.

Conclusions
The innovation of this paper was the joint use of Sentinel-2 and agrometeorological data in the application of the SAFER model, without thermal bands, which allowed quantification and spatial analyses of the crop coefficient and biomass production in watersheds with varied land uses located in southeastern Brazil at a 10 m resolution.This was done by first modeling the ratio of the actual (ET a ) and the reference (ET 0 ) evapotranspiration, analogous to Kc, at the satellite overpass time by the so-called SAFER algorithm, and then modeling the biomass production using the Monteith model, calibrated to use the crop coefficient found in the previous step of this methodology.This methodology captures both the characteristics of the phenological phases of the plants and the soil moisture conditions, allowing, in a resolution of 10 m, monitoring of the evapotranspiration and use of water by the crop in a precision scale.The surface albedo values, NDVI, net radiation, and surface temperature obtained in this work were corroborated by other studies done in southeastern Brazil and in other places, with characteristics that justify their comparison.

Figure 2 .
Figure 2. Basic land use of watersheds, from the day of year (DOY) 75/2017.

Figure 2 .
Figure 2. Basic land use of watersheds, from the day of year (DOY) 75/2017.

Figure 3 .
Figure 3. Flowchart for the large-scale radiation balances (net radiation, RN), crop coefficient (ETa/ET0), and biomass production (BIO) by applying the SAFER (Simple Algorithm for Evapotranspiration Retrieving) and Monteith method using Sentinel-2 images together with agrometeorological data (solar radiation, RG; average air temperature, TA; and reference evapotranspiration, ET0).

Figure 3 .
Figure 3. Flowchart for the large-scale radiation balances (net radiation, RN), crop coefficient (ET a /ET 0 ), and biomass production (BIO) by applying the SAFER (Simple Algorithm for Evapotranspiration Retrieving) and Monteith method using Sentinel-2 images together with agrometeorological data (solar radiation, R G ; average air temperature, T A ; and reference evapotranspiration, ET 0 ).

Figure 4 .
Figure 4. Rainfall (P), average temperature (TA), and solar radiation (RG) at a daily scale between 16 March 2017 and 2 April 2018 in watersheds located in the municipality of Águas de Santa Bárbara, the central-western part of São Paulo State, Brazil.The temporal location of the DOYs of the present research is indicated by arrows (DOY = day of year).

Figure 4 .
Figure 4. Rainfall (P), average temperature (T A ), and solar radiation (R G ) at a daily scale between 16 March 2017 and 2 April 2018 in watersheds located in the municipality of Águas de Santa Bárbara, the central-western part of São Paulo State, Brazil.The temporal location of the DOYs of the present research is indicated by arrows (DOY = day of year).

Figure 5 .
Figure 5. Surface albedo at a daily scale between 2017 and 2018 in watersheds with varied land use (DOY = day of year).

Figure 5 .
Figure 5. Surface albedo at a daily scale between 2017 and 2018 in watersheds with varied land use (DOY = day of year).

Figure 6 .
Figure 6.Net radiation (MJ) at a daily scale between 2017 and 2018 in watersheds with varied land use (DOY = day of year).

Figure 6 .
Figure 6.Net radiation (MJ) at a daily scale between 2017 and 2018 in watersheds with varied land use (DOY = day of year).

Figure 7 .
Figure 7. Surface temperature (K) at a daily scale between 2017 and 2018 in watersheds with varied land use.(DOY = day of year).

Figure 7 .
Figure 7. Surface temperature (K) at a daily scale between 2017 and 2018 in watersheds with varied land use.(DOY = day of year).

Figure 8 .
Figure 8. NDVI between 2017 and 2018 in watersheds with varied land use.(DOY = day of year).

Figure 8 .
Figure 8. NDVI between 2017 and 2018 in watersheds with varied land use.(DOY = day of year).

Figure 9 .
Figure 9. Crop coefficient between 2017 and 2018 in watersheds with varied land use (DOY = day of year).

Figure 9 .
Figure 9. Crop coefficient between 2017 and 2018 in watersheds with varied land use (DOY = day of year).
175, and 235 presented the highest ETr values because they are after intense precipitations and occur during the phenological stages of sugarcane and Eucalyptus with the highest water consumption.ETr characterizes the water status in the root zones[59].The most important reason for the highest ETr values were prior rainy days, which make the soil more moist, but it also depends on stomatal regulation and plant adaptation to water scarcity conditions [60].Zhou and Zhou [61] concluded that air humidity, air temperature, and the available energy were the most important variables for the Kc variations in a reed marsh located in Northeast China.The region in the north of the study area occupied by silviculture reaches the highest ETr values in DOY 175, before the beginning of the harvest season.The sugarcane areas located in the west of the study area have their highest ETr in DOY 175, before harvest, when ETr decreases in this region, as seen in DOY 235.Forest and pasture areas show ETr stability.

Figure 10 .
Figure 10.Variation of the mean values of the crop coefficient (ETr) with the accumulated degreedays (DDac) in the sugarcane crops of watersheds in the municipality of Águas de Santa Bárbara, the central-western part of São Paulo State, Brazil.A basal temperature of 18 °C was considered, following Doorenbos, Kassan (1979).

Figure 10 .
Figure 10.Variation of the mean values of the crop coefficient (ET r ) with the accumulated degree-days (DDac) in the sugarcane crops of watersheds in the municipality of Águas de Santa Bárbara, the central-western part of São Paulo State, Brazil.A basal temperature of 18 • C was considered, following Doorenbos, Kassan (1979).

Table 1 .
Mean values and standard deviations at a daily scale of surface albedo between 2017 and 2018 by land use and day of year (DOY).

Table 2 .
Mean values and standard deviations at a daily scale of surface albedo between 2017 and 2018 by land use and day of year (DOY).

Table 3 .
Mean values and standard deviation at a daily scale of surface albedo between 2017 and 2018 by land use and day of year (DOY).

Table 3 .
Mean values and standard deviation at a daily scale of surface albedo between 2017 and 2018 by land use and day of year (DOY).

Table 4 .
Mean values and standard deviation of NDVI between 2017 and 2018 by land use and day of year (DOY).

Table 4 .
Mean values and standard deviation of NDVI between 2017 and 2018 by land use and day of year (DOY).

Table 5 .
Mean values and standard deviation of the crop coefficient between 2017 and 2018 by land use and day of year (DOY).

Table 5 .
Mean values and standard deviation of the crop coefficient between 2017 and 2018 by land use and day of year (DOY).

Table 6 .
Mean values and standard deviation of the crop coefficient between 2017 and 2018 by land use and day of year (DOY).

Table 6 .
Mean values and standard deviation of the crop coefficient between 2017 and 2018 by land use and day of year (DOY).