Using Remote Sensing to Quantify the Joint Effects of Climate and Land Use/Land Cover Changes on the Caatinga Biome of Northeast Brazilian

  Citation: Jardim, A.M.d.R.F.; Araújo Júnior, G.d.N.; Silva, M.V.d.; Santos, A.d.; Silva, J.L.B.d.; Pandorfi, H.; Oliveira-Júnior, J.F.d.; Teixeira, A.H.d.C.; Teodoro, P.E.; de Lima, J.L.M.P.; et al. Using Remote Sensing to Quantify the Joint Effects of Climate and Land Use/Land Cover Changes on the Caatinga Biome of Northeast Brazilian. Remote Sens. 2022,14, 1911. https://doi.org/ 10.3390/rs14081911 Academic Editors: Baojie He, Ayyoob Sharifi, Chi Feng and Jun Yang Received: 3 March 2022 Accepted: 1 April 2022 Published: 15 April 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). remote sensing Article Using Remote Sensing to Quantify the Joint Effects of Climate and Land Use/Land Cover Changes on the Caatinga Biome of Northeast Brazilian Alexandre Maniçoba da Rosa Ferraz Jardim 1,2,*, George do Nascimento Araújo Júnior 1,2, Marcos Vinícius da Silva 1, Anderson dos Santos 1, Jhon Lennon Bezerra da Silva 1, Héliton Pandorfi 1, JoséFrancisco de Oliveira-Júnior 3, Antônio Heriberto de Castro Teixeira 4, Paulo Eduardo Teodoro 5, João L. M. P. de Lima 6,7 , Carlos Antonio da Silva Junior 8, Luciana Sandra Bastos de Souza 2, Emanuel Araújo Silva 9and Thieres George Freire da Silva 1,2 1Department of Agricultural Engineering, Federal Rural University of Pernambuco, Recife 52171-900, Brazil; [email protected] (G.d.N.A.J.); mar[email protected] (M.V.d.S.); [email protected] (A.d.S.); [email protected] (J.L.B.d.S.); [email protected] (H.P.); [email protected] (T.G.F.d.S.) 2Academic Unit of Serra Talhada, Federal Rural University of Pernambuco, Serra Talhada 56909-535, Brazil; [email protected] 3Institute of Atmospheric Sciences, Federal University of Alagoas, Maceió57072-970, Brazil; [email protected] 4Water Resources Department, Federal University of Sergipe, São Cristóvão 49100-000, Brazil; [email protected] 5Department of Agronomy, Federal University of Mato Grosso do Sul, Chapadão do Sul 79560-000, Brazil; [email protected] 6MARE—Marine and Environmental Sciences Centre, University of Coimbra, 3000-456 Coimbra, Portugal; [email protected] 7Department of Civil Engineering, Faculty of Sciences and Technology, University of Coimbra, 3030-788 Coimbra, Portugal 8Department of Geography, State University of Mato Grosso (UNEMAT), Sinop 78555-000, Brazil; [email protected] 9Department of Forest Sciences, Federal Rural University of Pernambuco, Recife 52171-900, Brazil; [email protected] *Correspondence: alexandre.jar[email protected] Abstract: Caatinga biome, located in the Brazilian semi-arid region, is the most populous semi-arid region in the world, causing intensification in land degradation and loss of biodiversity over time. The main objective of this paper is to determine and analyze the changes in land cover and use, over time, on the biophysical parameters in the Caatinga biome in the semi-arid region of Brazil using remote sensing. Landsat-8 images were used, along with the Surface Energy Balance Algorithm for Land (SEBAL) in the Google Earth Engine platform, from 2013 to 2019, through spatiotemporal modeling of vegetation indices, i.e., leaf area index (LAI) and vegetation cover (V C ). Moreover, land surface temperature (LST) and actual evapotranspiration (ET a ) in Petrolina, the semi-arid region of Brazil, was used. The principal component analysis was used to select descriptive variables and multiple regression analysis to predict ET a . The results indicated significant effects of land use and land cover changes on energy balances over time. In 2013, 70.2% of the study area was composed of Caatinga, while the lowest percentages were identified in 2015 (67.8%) and 2017 (68.7%). Rainfall records in 2013 ranged from 270 to 480 mm, with values higher than 410 mm in 46.5% of the study area, concentrated in the northern part of the municipality. On the other hand, in 2017 the lowest annual rainfall values (from 200 to 340 mm) occurred. Low vegetation cover rate was observed by LAI and V C values, with a range of 0 to 25% vegetation cover in 52.3% of the area, which exposes the effects of the dry season on vegetation. The highest LST was mainly found in urban areas and/or exposed soil. In 2013, 40.5% of the region’s area had LST between 48.0 and 52.0 ◦ C, raising ET a rates (~4.7 mm day −1 ). Our model has shown good outcomes in terms of accuracy and concordance (coefficient of determination = 0.98, root mean square error = 0.498, and Lin’s concordance correlation Remote Sens. 2022,14, 1911. https://doi.org/10.3390/rs14081911 https://www.mdpi.com/journal/remotesensing Remote Sens. 2022,14, 1911 2 of 27 coefficient = 0.907). The significant increase in agricultural areas has resulted in the progressive reduction of the Caatinga biome. Therefore, mitigation and sustainable planning is vital to decrease the impacts of anthropic actions. Keywords: tropical dry forest; surface energy balance; Brazilian semi-arid; SEBAL; actual evapotranspiration 1. Introduction The Caatinga biome occupies a large portion of the Brazilian semi-arid region. It has a high ecological diversity of plant species, and it is considered the largest in the world under semi-arid conditions [ 1 , 2 ]. It occupies an area of 900,000 km 2 , which corresponds to approximately 70% of the northeast region of Brazil (NEB). However, only 7.5% of this habitat is protected by law [ 3 – 5 ]. The municipality of Petrolina, PE, Brazil, is located in the semi-arid region of the Caatinga biome. It is within the hydrographic basin of the São Francisco river, which favors the development of irrigated agriculture in the region, mainly fruit-growing (e.g., grapes and mangoes) [ 6 – 8 ], and it leads to positive socioeconomic implications, but also increases conflicts regarding water use [2,9]. Caatinga is located in the world’s most populated dry area, with more than 53 million inhabitants and a population density close to 34 inhabitants per km 2 . This biome has been affected by environmental degradation and loss of biodiversity over the last decades, mainly by the intensification of agriculture (e.g., rainfed and irrigated crops cultivation), urban expansion, and the advance of pasture areas replacing the natural vegetation [ 10 – 12 ]. Expanding agricultural activities has led to the deforestation of native areas, soil disturbance, changes in the hydrological cycle, and higher carbon emissions [ 13 – 15 ]. Together with climate changes, such as reduced rainfall and intensified drought events, this makes the Caatinga biome and the Brazilian ecosystem the most threatened and susceptible to desertification. Furthermore, these changes have been increasingly compromising natural resources and environmental sustainability [ 16 – 19 ] due to the changes in surface properties and biophysical variables, such as vegetation cover (V C ), land surface temperature (LST), and evapotranspiration (ET) [ 20 – 22 ]. ET is one of the main response parameters of vegetated areas as a function of local water conditions, and it is also an important component of the hydrological cycle [23,24]. Therefore, monitoring physical–water indicators of environmental change conditions, such as the loss of biodiversity of biomes and land use and occupation, is vital in managing scarce water resources. Furthermore, these indicators may be helpful in the planning of agricultural activities, as well as in the management of dry areas and the sustainable use and management of natural resources [17,25–30]. In this scenario, remote sensing has been used as a tool that presents fast and low operational costs, being efficient in calculating the biophysical parameters used in the energy, water, and vegetation balances (e.g., [ 31 – 33 ]). In recent years, remote sensing has also been considered a good alternative for replacing expensive and difficult-to-obtain equipment used in in situ studies [ 34 ]. Furthermore, the modeling used in these balances is performed with the help of algorithms, which are essential techniques in the extraction of information from satellite images on a regional and global scale [ 27 , 28 , 33 , 35 ]. It may have its efficiency improved by the use of open-source cloud programming languages, e.g., Google Earth Engine (GEE). In GEE, the user can effectively implement algorithms and process large data volumes [36]. In this context, there are several surface energy balance models, such as Mapping EvapoTranspiration at high Resolution with Internalized Calibration (METRIC), Surface Energy Balance System (SEBS), Simple Algorithm for Evapotranspiration Retrieving (SAFER), Atmosphere–Land Exchange Inverse (Alexi), and the Energy Balance Algorithm for Land (SEBAL), based on remote sensing data. The approach is substantiated on biophysical parameters, such as LST, albedo, Normalized Difference Vegetation Index (NDVI), and Remote Sens. 2022,14, 1911 3 of 27 emissivity, that are crucial for estimating ET [ 16 , 37 – 42 ]. The SEBAL algorithm has been widely and successfully applied to various world ecosystems, including semi-arid conditions in Brazil [ 43 – 46 ]. Moreover, this method seeks to eliminate the propagation of errors in the partitioning of the energy balance and the need for atmospheric correction in the estimate of surface temperature. These interactions allow the generation of the sensible heat flux corrected for atmospheric stability and instability conditions [ 47 – 51 ]. NDVI is one of the most widely used indices in the literature for vegetation cover analysis, environmental degradation, and vegetation primary production resilience, and it also helps in monitoring agricultural crops, forest ecosystems, and drought assessments [ 52 – 55 ]. With applications in several countries, Bastiaanssen et al. [ 49 ] have obtained highly accurate results using the SEBAL algorithm on different vegetated surfaces with forests, agricultural crops (irrigated and rainfed), and even extreme landscapes, such as desert areas. Previous studies have shown that the SEBAL can be used in Petrolina, PE, Brazil, for the algorithm has already been calibrated and validated for the region in various ecosystems with good agreement between orbital images and field measurements [44,56–59]. Based on the above, the main objective of this paper is to determine and analyze the changes in land cover, land use, and occupation on biophysical parameters in the Caatinga biome in the semi-arid region of Brazil. In the present study, the changes were assessed from 2013 to 2019 using Landsat imagery, and biophysical parameters (net radiation, energy balance, LAI, V C , ET a , and LST) were estimated by utilizing remote sensing. Additionally, the SEBAL model was applied to determine the surface energy balance in different vegetation environments. After applying the model, the dataset provided by SEBAL allowed us to determine the turbulent fluxes and actual evapotranspiration (ET a ) of the different land use and land cover (LULC) types. Thus, through the SEBAL products, the spatiotemporal variation patterns of ETain agricultural and forestry areas were evaluated. 2. Materials and Methods 2.1. Study Area The study was carried out in the municipality of Petrolina, located in the State of Pernambuco, Brazil. The region comprises the domain of the Caatinga biome and it belongs to the semi-arid region of the sub-mean of the São Francisco Valley. The town is considered the largest fruit-growing center of the Brazilian semi-arid region due to the easy access to the São Francisco river, which supplies the irrigated perimeters. The municipality comprises a territorial area of 4561.870 km 2 (Figure 1), with an estimated population of 349,145 inhabitants [38]. A characteristic of the Caatinga biome is the presence of different floristic mosaics, consisting of an area of tree and shrub vegetation, which presents its distribution conditioned to climatic and environmental variations, especially rainfall intensity and frequency [ 60 ], as well as geological configurations and soil properties [ 5 ]. The vegetation of this biome presents deciduous species adapted to water deficit conditions and with expressive biomass production in rainy seasons, resulting from the local climatic conditions [ 40 , 61 ]. The canopy cover of the species of the Caatinga biome presents discontinuous characteristics, making possible the soil exposure in dry periods, presence of herbaceous stratum, cactus species, and shrubs [ 23 , 61 ]. It is worth noting that this region presents an expressive modification of the native landscape (Caatinga) in areas of irrigated agricultural cultivation. Remote Sens. 2022,14, 1911 4 of 27 Figure 1. Spatial location of the study area, municipality of Petrolina, Pernambuco, Northeast Brazil. According to the Köppen–Geiger climate classification, the region’s climate is of the BSh type, characterized as semi-arid tropical, with an average air temperature of 26.4 ◦ C, average relative humidity of 62%, and annual rainfall of 520 mm [ 62 , 63 ]. Rainfall pattern is irregular throughout the year, resulting from its geographical location and the Intertropical Convergence Zone (ITCZ) influence, with rainfall predominating from February to May [ 10 , 64 , 65 ]. The predominant soils in the municipality are Typic Quartzipsamment, Ultisol Plinthic, Arenosol, and Haplic Acrisol [66–68]. 2.2. Satellite Images and Weather Datasets Annual records of average air temperature ( ◦ C), global radiation (MJ m −2 ), relative humidity (%), atmospheric pressure (kPa), wind speed (m s −1 ), and rainfall (mm) were obtained from the database of the National Institute of Meteorology [ 69 ] (Figure 2). The average air temperature (T a ) ranged between 24.1 and 30.7 ◦ C, in which 2015 and 2019 were the warmest years studied (T a 28 ◦ C). Overall, November, December, January, February, and March presented Tavalues above 28 ◦C (Figure 2). The years 2013 (334.4 mm) and 2019 (221.6 mm) showed the highest rainfall rates. Most of the rainfall recorded for Petrolina was concentrated from December to March (Figure 2). Such climatic conditions concerning the municipality were also reported in the literature (e.g., [2,70]). Rainfall data corresponding to the 30 days prior to the four imaging dates studied used were obtained from the Climate Hazards Group InfraRed Precipitation with Station (CHIRPS). CHIRPS are new precipitation products covering the coordinates 50 ◦ S–50 ◦ N and 180 ◦ E–180 ◦ W, with 0.05 ◦ ( ± 5.3 km) spatial resolution and daily to seasonal, temporal resolutions, available worldwide since 1981 [ 71 ]. CHIRPS data were extracted from the Google Earth Engine platform (https://earthengine.google.com/, accessed on 20 August 2021) using JavaScript programming language. Then, they were exported in spreadsheet format (*.xls), using the dataset since 1981 from the collection ee.ImageCollection (“UCSBCHG/CHIRPS/DAILY”). Remote Sens. 2022,14, 1911 5 of 27 Figure 2. Monthly meteorological variations (rainfall, average air temperature, and global solar radiation) of the municipality of Petrolina, Brazil, from 2013 to 2019. We used Operational Land Imager (OLI) Collection 1 Level 1 bands 2 (0.450–0.51 µ m), 3 (0.53–0.59 µ m), and 4 (0.64–0.67 µ m) in the visible spectrum, 5 (0.85–0.88 µ m) in the near-infrared, and 6 (1.57–1.65 µ m) and 7 (2.11–2.29 µ m) in the shortwave infrared, all with a spatial resolution of 30 m, as well as band 10 from the Thermal Infrared Sensor (TIRS) with a 100 m spatial resolution. Besides this, we used four Landsat-8 OLI/TIRS images, path 217 and row 66, corresponding to the years 2013, 2015, 2017, and 2019 (for the dates and times of the satellite overpass, see Table 1). The choice criteria adopted were the absence of clouds (10%) and the images corresponding to the transition period between the dry and rainy seasons in the region under study, from the years of 2013 to 2019. This period presented severe and extreme drought events in the northeast [ 23 , 72 ]. All the images were obtained from the United States Geological Survey (USGS) platform (https://earthexplorer.usgs.gov/, accessed on 10 August 2021) and processed through the Land Surface Reflectance Code (LaSRC). Table 1. Date of the Landsat-8 satellite pass, followed by the Julian day (JD), Earth–Sun distance (dr, astronomical units—AU), local time of the equator pass (h, hour; min, minutes), zenith angle ( θ , ◦ ), solar elevation angle (E, ◦ ), and sun azimuth angle ( ϕ , ◦ ) for the municipality of Petrolina, Pernambuco, Brazil. Acquisition Date JD dr Local Time θEϕ 5 October 2013 278 0.99 9 h 49 min a.m. 0.90 65.12 82.93 12 November 2015 316 0.99 9 h 48 min a.m. 0.90 64.85 113.46 16 October 2017 289 0.99 9 h 48 min a.m. 0.91 65.81 92.82 7 November 2019 311 0.99 9 h 48 min a.m. 0.90 65.41 110.36 Note: zenith angle (θ) = sin(E). Source: USGS/NASA [73]. 2.3. Vegetation Indices The Normalized Difference Vegetation Index (NDVI) was calculated to represent the amount and quality of vegetation present on the surface, characterized as an indicator of wet conditions, calculated using Equation (1). NDVI =ρNIR −ρRed ρNIR+ρRed (1) where ρNIR and ρRed are the reflectances measured in the near-infrared and red bands (i.e., Landsat-8 multispectral bands 5 and 4 of the OLI sensor), respectively, they range from − 1 to +1. Values close to 1 on a positive scale correspond to high photosynthetic activity, and when negative, generally correspond to water bodies. Remote Sens. 2022,14, 1911 6 of 27 Based on the NDVI, we calculated the vegetation cover (V C ) of the study area (Equation (2)) , according to Gao et al. [54]. VC=NDVI −NDVIS NDVIV−NDVIS ·100 (2) where V C is the vegetation cover, NDVI S is the minimum NDVI value from bare soil pixels obtained in the study area, and NDVI V is the maximum NDVI value found in vegetated areas, i.e., from fully vegetated pixels. The NDVI S and NDVI V used to calculate V C were obtained from the domain of each NDVI image percentile map obtained from the NDVI histograms. Soil-Adjusted Vegetation Index (SAVI) was calculated to observe the vegetation cover of the area (Equation (3)). SAVI =(1+L)·(ρNIR −ρRed) (L+ρNIR+ρRed)(3) where Lis the adjustment factor to the soil, which varies between 0 and 1. The value 0 does not reach change, and resembles the NDVI. In areas with low-density vegetation, the value 1 is assigned; for intermediate-density vegetation areas, the value of 0.5; and for areas with high-density vegetation, the value 0.25 is assigned [ 74 ]. The adjustment factor of 0.5 was adopted due to the study region indicating an intermediate vegetation coverage in most of the year, with predominant vegetation of the Caatinga biome, in the Brazilian semi-arid region [75–77]. To evaluate changes in vegetation biomass, the leaf area index (LAI, m 2 m −2 ) was determined (Equation (4)), a fundamental biophysical variable for monitoring studies of agricultural land and vegetation moisture conditions [50]. LAI = −ln0.69−SAVI 0.59  0.91 (4) 2.4. Methodology for Estimating Evapotranspiration Using Satellite Images Evapotranspiration was estimated using the Energy Balance Algorithm for Land (SEBAL). For this, routine meteorological data and spectral bands from the Landsat-8 satellite were used. The SEBAL algorithm was implemented using JavaScript code through the Google Earth Engine (GEE) platform. The data were exported in spreadsheet format (*.xls). SEBAL uses mathematical modeling and operations to calculate the surface energy balance components and determine evapotranspiration. Thus, energy balance components are computed pixel-by-pixel, as described in Figure 3. Here, we show that the algorithm has a good performance and high accuracy, as well as being calibrated and validated with simultaneous field and Landsat satellite measurements [47–49,56,59,61,78,79]. Remote Sens. 2022,14, 1911 7 of 27 Figure 3. Flowchart of the SEBAL model for estimating evapotranspiration. Note: LST and LSE are the land surface temperature and land surface emissivity, respectively, NDVI is the Normalized Difference Vegetation Index, εa is the atmospheric emissivity, u* is the friction velocity, T a is the air temperature, R n is the net radiation, Gis the soil heat flux, z om is the momentum roughness length, z 1 and z 2 are the two heights between the surface of the anchor pixels, r ah is the near-surface aerodynamic resistance to heat transport, DEM is the digital elevation model, dT is the near-surface air temperature gradient, aand bare the calibration coefficients, LST or T s is the land surface temperature, ρair is the air density, C p is the specific heat of air, His the sensible heat flux, ψm and ψh are the stability correction factors for momentum and sensible heat, respectively, kis the von Karman constant, u 200 is the wind speed at the height of 200 m, and Λ is the evaporative fraction. These preand processing steps were performed inside Google Earth Engine (GEE) cloud platform. The acronyms and symbols used in this study are summarized in the Abbreviations section. Surface Albedo Adjustment The surface albedo ( αsup ) corresponds to a measure of the reflectivity of the Earth’s surface, for each pixel, with atmospheric correction obtained according to Equation (5) [ 50 , 80 ]. αsup =αtoa −αpath τsw2(5) where αtoa is the albedo at the top of the atmosphere, that is, before atmospheric correction, αpath is the atmospheric reflectance (set to 0.03, as used by Silva et al. [ 80 ]), and τsw is the atmospheric transmissivity for clear sky conditions, according to Equation (6) [50,80]: τsw=0.35 +0.627 ·exp"−0.00146 ·Pa Kt·cos(θ)−0.075W cos(θ)0.4#(6) where P a is the atmospheric pressure (kPa), with dataset available freely in https://portal. inmet.gov.br/ (accessed on 10 August 2021), K t is the turbidity coefficient of the atmosphere (K t = 1.0, for a clear sky day), according to Allen et al. [ 47 ] and Silva et al. [ 80 ], θ is the solar zenith angle, and Wis the precipitable water (mm), estimated from Equation (7) [81]. W=0.14 ·ea·Pa+2.1 (7) where e a is the actual atmospheric water vapor pressure (kPa), estimated from Equation (8). ea=HR ·es 100 (8) Remote Sens. 2022,14, 1911 8 of 27 where HR is the instantaneous relative humidity (%), and e s is the water vapor saturation pressure (kPa), estimated from Equation (9). es=0.6108 ·exp17.27 ·T0 237.3 +T0(9) where T0is the instantaneous air temperature (◦C) at the moment of the satellite pass. For obtaining αtoa , a linear combination of the spectral reflectance of the six reflective OLI bands was performed according to Equation (10) [80]: αtoa =0.300 r2+0.277 r3+0.233 r4+0.143 r5+0.036 r6+0.001 r7(10) where r 2 , r 3, r 4 , r 5, r 6 , and r 7 are the surface spectral reflectances for bands 2, 3, 4, 5, 6, and 7 of the Landsat-8 OLI, respectively. We use Equation (11) to obtain each of the spectral reflectances. rb=Addb+Multb·DN cos(θ)·dr (11) where the terms Add b and Mult b belong to the radiometric rescaling group, specifically reflectance_add_band (equal to − 0.1) and reflectance_mult_band (equal to 0.00002), respectively, presented in the metadata of each OLI—Landsat-8 image, DN is the digital number value corresponding to the pixel, θ is the solar zenith angle at the data acquisition time, and dr is the Earth–Sun distance in astronomical units. 2.5. Determination of Surface-energy Partitioning Based on the surface energy balance components, the evaporative fraction was determined. Initially, the surface radiation balance or net radiation—R n was calculated, which is distributed by the energy partitioning in front of the sensible heat fluxes—H, latent—LE, and soil heat flux—G[ 47 – 50 , 61 , 78 , 79 ]. By performing this process, a linear relationship between the surface and air temperature gradient was considered to exist. From this relationship and the internal calibration process for extreme conditions such as temperature and humidity, it was established the need for obtaining the knowledge of the so-called “anchor pixels”, i.e., hot and cold pixels, which are indicative of zero and maximum evapotranspiration, respectively [47,48,61] (Figure 3). The land surface temperature (LST) in K (Kelvin) was obtained using the spectral radiance in band 10 of the TIRS sensor and the emissivity in the nearest band— εnb by the modified Planck’s Law [82], as described in Equation (12). LST =K2 ln” nb·K1 L10 +1(12) where K 1 and K 2 are radiation constants specific for the Landsat-8 TIRS band 10, equaling 774.89 W m−2sr−1µm−1and 1321.08 K, respectively, provided by NASA/USGS; and L10 is the radiance at the wavelength received by the sensors (band 10, the thermal band). The εnb was calculated based on the LAI for each pixel according to Equation (13) [ 41 ]. ” nb=0.97 +0.0033 ·LAI (13) Initially, the temperature variation and aerodynamic resistance to heat transport in all pixels of the study area (Petrolina, Pernambuco) were determined. The atmosphere was initially assumed to be in a neutral stability condition. For this study, the hot pixel was considered in the exposed soil plots (i.e., no vegetation cover and/or little vegetation and low moisture content), assuming LE equal to zero. The cold pixel was considered in grape orchard plots irrigated by micro-sprinklers, when Hcan be considered zero [ 47 , 48 , 61 , 78 , 79 , 83 ] (see Figure 3). Since turbulent effects affect atmospheric conditions and air resistance, the Remote Sens. 2022,14, 1911 9 of 27 Monin–Obukhov similarity theory was applied and considered in the computation of Hin all pixels of the study area. It is worth noting that the Monin–Obukhov length was used for corrections to the initial stable condition of the atmosphere [47–49,78,79]. Calculation of Energy Fluxes (Hand LE) and Evaporative Fraction The sensible heat flux (H) in SEBAL is calculated using an iterative procedure from the aerodynamic function (Equation (14)) [47–49,78,79]. H=ρair ·Cp·(a+b·LST) rah (14) where ρair is the moist air density (kg m −3 ), C p is the air specific heat at constant pressure (1004 J kg −1 K −1 ), aand bare calibration constants of the temperature difference between two heights (i.e., between the roughness length for heat transfer and the reference height, usually 0.1 and 2.0 m above the displacement plane), and r ah is the near-surface aerodynamic resistance to heat transport (s m −1 ). Fundamentally, the coefficients aand bare determined through an internal calibration for each satellite image by interactive processes. We consider extreme pixels of wet/cold and dry/hot spots. They were selected to develop a linear relationship between the aerodynamic temperature of the surface and the air temperature difference, and the LST. By knowing the components of the surface energy balance, such as the net radiation (R n ,Wm −2 ), sensible heat flux (H,Wm −2 ), and soil heat flux (G,Wm −2 ), the latent heat flux (LE,Wm −2 ) was determined, both corresponding to the time of the satellite pass over the study area, according to Equation (15). LE =Rn−H−G(15) Subsequently, we determined the evaporative fraction ( Λ ) according to Equation (16). Λ=LE Rn−G(16) 2.6. Estimate of ETaUsing SEBAL Method Finally, as a SEBAL product, we determine the actual evapotranspiration (ET a , mm day −1 ) [ 84 ] based on Equation (17) below, for each satellite image used. In this study, the implemented SEBAL model had already been extensively validated and calibrated under forests and agricultural land conditions [56,61,85]. ETa=Λ·Rn24 ·86, 400 ˘(17) where R n24 is the daily net radiation (W m −2 ), 86,400 is a constant for daily timescale conversion (i.e., converts from seconds to days), and λ is the latent heat of vaporization of water (J kg −1 ). Then, the latent heat of vaporization allows the ET a expression in mm day −1 . Hence, accurate estimation of R n24 (Equation (18)) was determined according to Bastiaanssen et al. [49], and Lee and Kim [86]: Rn24=Λ·(1−αsup)·Rn−a·τsw(18) where ais a regression coefficient of the relationship between net longwave radiation and atmospheric transmissivity on a daily scale, to which we assigned the value 143, as proposed by Teixeira et al. [ 61 ]. The acronyms and symbols used in this study are summarized in the Abbreviations section. Remote Sens. 2022,14, 1911 16 of 27 distribution of rainfall accumulation in the previous days. On the other hand, during the dates studied here, the highest mean LAI values stood out only in areas of arboreal Caatinga (0.58 ± 0.45 m 2 m −2 ) and agriculture (0.75 ± 0.49 m 2 m −2 ), with the maximum values being associated with irrigated agricultural areas, more specifically orchards, with increased biomass production. However, in areas of pasture, urban infrastructure, and areas of arboreal and herbaceous Caatinga, the LAI values were close to or equal to zero, making its variation more homogeneous, that is, closer to the daily mean, with standard deviation (SD) ranging between 0.05 and 0.10 m2m−2within the land cover classes (Figure 7). 3.4. Land Surface Temperature (LST) in the Studied Classes In the present study, the minimum LST values were seen in areas of water bodies, while the maximum LST values were seen in areas of exposed soils, located at points of pasture and degraded Caatinga, urban infrastructure (asphalt, concrete, and gravel surfaces), and agricultural areas undergoing soil preparation for cultivation (Figure 8). This result is expected in bare soil locations under intense anthropic activity due to the transformation of land use/land cover classes into non-evaporating surfaces. This makes the place’s temperature higher and reduces water availability in the soil, which causes serious problems in agricultural crops. The computed LST map is shown in Figure 8. Figure 8. Spatiotemporal distribution of the land surface temperature—LST ( ◦ C) in the municipality of Petrolina, Pernambuco, Brazil, on the imaging dates 5 October 2013 ( a ), 12 November 2015 ( b ), 16 October 2017 (c) and 7 November 2019 (d). Due to the replacement of primary vegetation with pastures, agricultural crops, and urban occupation, changes in land use can substantially affect the heat and mass exchange in the soil–plant–atmosphere system, propitiating the retention of a higher amount of heat by the Earth’s surface [ 2 , 52 ]. Land abandonment and excessive mechanical disturbance of the soil may also alter the heat exchange with the environment and cause lower thermal and radiant energy lag; thus, the land conversion had increased LST in the area of the non-evaporating surfaces. It can be observed that the average LST values for the dates studied were higher in the areas dominated by pasture (47.69 ± 1.47 ◦ C), herbaceous Caatinga (47.28 ± 1.27 ◦ C), and shrub Caatinga (46.07 ± 1.44 ◦ C), even higher than those observed in areas with urban Remote Sens. 2022,14, 1911 17 of 27 infrastructure (45.80 ± 1.52 ◦ C) (Figure 8). According to Zhao et al. [ 110 , 111 ], sites with different land cover types may have an LST increase gradient along the urban to rural profile. The high LST in the pasture, herbaceous, and shrub Caatinga areas is related to the lower percentage of ground covering by vegetation (see Figure 6), which results in drier exposed soil, with higher albedos and lower evaporative cooling flux rates, a factor that increases LST. Another related factor contributing to the high LST of pastures is that grasses have shallower roots. Therefore, they can only access the water available in the superficial soil layers, which depletes faster than in deeper layers [ 2 ]. Vegetation canopy can retain rainwater and decrease groundwater recharge by altering evaporative flux and raising the land surface temperature. On the other hand, it is observed, in general, that in areas dominated by arboreal Caatinga and agriculture, the average LST values are lower (38.99 ± 2.48 ◦ C and 43.11 ± 2.30 ◦ C, respectively) (Figure 8). The highest vegetation cover and the highest soil humidity in the areas dominated by these classes favor LST reduction. However, they present the greatest spatial variations of LST among all land use and land cover classes, according to the standard deviation values ( ± SD). The main advantage of using LST data from satellite images is the total surface coverage. In this way, each time series of pixels of the LST map can be considered a “virtual weather station” [112]. 3.5. Variations of the Actual Evapotranspiration (ETa) of Land Use Classes Rainfall regime directly influenced ET a , so that on 16 October 2017 (Figure 9c), the date with the lowest rainfall accumulation in the previous days, the lowest ET a rates were found, with an average value of 2.02 mm day −1 . On the other hand, on 7 November 2019 (Figure 9d), the period with the highest rainfall accumulation, average ET a rates were 2.62 mm day −1 (Figure 9). According to Teixeira et al. [ 20 ], high evapotranspiration values in Caatinga areas occur right after rains. Therefore, the previous rainfall raises soil water availability and keeps native species with turgid structures and greener canopy. Figure 9. Spatiotemporal distribution of actual evapotranspiration (ET a , mm day −1 ), calculated with Surface Energy Balance Algorithm for Land (SEBAL), in the municipality of Petrolina, Pernambuco, Brazil, on the imaging dates 5 October 2013 ( a ), 12 November 2015 ( b ), 16 October 2017 ( c ), and 7 November 2019 (d). Remote Sens. 2022,14, 1911 18 of 27 The ET a estimated by the SEBAL model showed variation both within and between land use and land cover classes. The lowest ET a observations in all the evaluated dates were observed in the areas occupied by pasture and mosaic of agriculture and pasture classes, with average values of 0.70 ± 0.73 mm day −1 and 1.01 ± 0.96 mm day −1 , respectively (Figure 9). These areas have dryland cultivation practices, which causes the lower water availability to affect the evapotranspiration rates; moreover, the heterogeneity of the areas causes sudden variations in ETa. On the other hand, due to the effect of the increase in air temperature and atmospheric demand verified throughout the dry season in Petrolina, combined with the presence of preserved riparian forests along stretches of water bodies and continuous irrigation in crops, higher mean ET a values were observed in the arboreal Caatinga (4.73 ± 0.49 mm day −1 ) and agriculture (3.07 ± 1.23 mm day −1 ) classes. For presenting areas with irrigated and dry cultivation, the average values of this class become more variable. When there is a greater contribution of moisture added to more dense vegetation, there is a favoring of the local microclimate in the region [113,114], a phenomenon reported in areas of arboreal vegetation. The shrub and herbaceous Caatinga classes showed greater heterogeneity indicated by the largest standard deviations (2.42 ± 0.76 mm day −1 ) and (1.46 ± 0.71 mm day −1 ), respectively, relative to the arboreal Caatinga class (Figure 9). Folhes et al. [ 115 ] reported that the evapotranspiration values (2.0 mm day −1 ) during the dry season for species of herbaceous–shrubby Caatinga. In dry periods, the Caatinga vegetation uses the available energy as sensible heat flux (H), limiting transpiration and photosynthesis, thus reducing evapotranspiration values [ 20 ]. However, arboreal vegetation is able to compensate the high vapor pressure deficit in the air, even in dry periods, when compared to shrub and herbaceous vegetation, due to the deep root system keeping up with the water stored in the soil [70,116,117]. The arboreal and shrub species play a fundamental eco-hydrological role, maintaining soil humidity and structuring its porosity, guaranteeing the maintenance of infiltration capacity and favoring the survival of species [116]. 3.6. Statistical Relations between the Variables Studied and Land Use In this study, we performed a PCA of the environmental variables in relation to land use and land cover classes. Therefore, the first two components with eigenvalues greater than 1.0 were extracted separately for the years 2013, 2015, 2017, and 2019 (Figure 10). In 2013, the two principal components explained 94.77% of the total variation, with 70.39% in the principal component 1 (PC1) and 24.38% in the principal component 2 (PC2). On the other hand, in 2015, 2017, and 2019, when added together, PC1 and PC2 represented 94.01, 92.34, and 94.56% of the total variation, respectively (Figure 10). In addition, it can be seen that the LULC class, with the least influence on components 1 and 2 in all years, is shrub Caatinga, with average eigenvalues (0.31 and 0.37, respectively). In addition, the classes with the greatest influence on PC1 with positive and negative eigenvalues are water bodies (2.69), agriculture ( − 2.35), arboreal Caatinga ( − 1.92), pasture (0.34), mosaic (0.16), and urban area (0.35). In PC2, the classes of LULC were water bodies ( − 4.36), mosaic (2.10), agriculture (−1.0), and arboreal Caatinga (−3.40) (Figure 10). Remote Sens. 2022,14, 1911 19 of 27 Figure 10. Scores obtained by principal component analysis (PCA) of environmental variables and land use and land cover. PC1 and PC2 are the first and second dimensions of PCA data, respectively. The four inserted panels below the PCA scores plots refer to the loadings plots of the first two principal components from 2013 to 2019. Through PCA, we observed that the ordering of variables in each principal component (PC) of the axes was influenced by the degree of vegetation cover and surface water status of the LULC classes (Figure 10). Thus, PC1 contributed more to the variability of the response of variables related to energy balance. Therefore, in PC2, the land use and land cover classes (i.e., pasture, mosaic, urban area, and herbaceous Caatinga) influenced the variables R n ,LE, ET a , and emissivity with higher mean loadings ( − 0.95, − 0.94, − 0.99, and − 0.80, respectively). On the other hand, they showed a high correlation with the variables H, LST, and albedo. For these LULC classes, this may be related to the presence of bare soil and thin vegetation covers that affect the regional microclimate and soil–plant–atmosphere system fluxes, resulting in higher albedos, lower rates of evaporative cooling fluxes, and highest LST [2,105]. On the other hand, PC2 contributed more to the variability of the responses of the variables regarding the canopy interactions (LAI and V C ) due to the strong negative correlation with cultivated land (i.e., agriculture class) (Figure 10). In all years, there was a predominance for LE and LAI in areas with arboreal Caatinga and agriculture, with agriculture presenting the highest V C . Caatinga presents a strong relationship with LAI in rainy periods due to the greater availability of water in the soil [ 118 ]. In the present study, the samples were taken in the period with low rainfall, so the LAI was not expressive compared to agricultural areas, which use irrigation and therefore increase the LAI. Thus, Remote Sens. 2022,14, 1911 20 of 27 the emissivity in the Caatinga vegetation is lower, for the emissivity of the soil is generally lower than that of the leaves [2]. The results of the analysis of variance (ANOVA) and multiple regression analysis of the established model are presented in Table 3. For the combination of two-variable models, the two best pre-established parameters based on the PCA results were the variables LST and H, used to build the regression model. In particular, these two variables are of great relevance in transferring energy to the atmosphere. The ANOVA results also showed that the model values are significant. Due to the observations, a joint analysis of the coefficient of determination obtained (R 2 = 0.98) can be performed, emphasizing the P-value obtained from our regression model, which was less than 0.001, thus indicating greater model accuracy and reliability (Table 3). Notably, it can be seen that the model’s F-value was 16,692.84, being greater than the critical value of F 0.05 = 3.018, which confirms the significance of the proposed model. In addition, the variables used provided a high R 2 and LCCC, being essential for the model’s accuracy. Furthermore, the results showed that the multiple linear regression model to determine ET a achieved a coefficient of determination of 0.98, RMSE of 0.498, MAE of 0.413, and dequal to 0.9620. It also resulted in PBIAS, NSE, and LCCC values averaged between − 13.32%, 0.826, and 0.907, respectively (Table 3). From the statistical analysis, this model is described as ET a = 6.89 − 0.0527LST − 0.0120H. Based on the RMSE, LCCC, and dof this application, the use of the ET a model, besides presenting a strong correlation between the variables Hand LST, as seen in Figure 10, expresses a biophysical model with enough efficiency and high agreement to determine ETa. Table 3. Analysis of variance (ANOVA) and regression coefficients results for the suggested model. Source of Variation df SS MS F-Value p-Value Regression 2 464.41 232.21 16,692.84 0.0001 LST 1 335.29 335.29 24,103.8 0.0001 H1 129.12 129.12 9281.9 0.0001 Error 397 5.52 0.01 Total 399 469.93 Regression statistics Predictors in model Regression coefficients β0β1β2R2 LST, H6.89 −0.0527 −0.0120 0.98 Model Statistical metrics RMSE MAE PBIAS (%) NSE LCCC d 0.498 0.413 −13.32 0.826 0.907 0.9620 df: degrees of freedom, SS: sum of squares, MS: mean square, LST: land surface temperature, H: sensible heat flux, β0 : intercept, β1 and β2 : estimated coefficient for the factor x, R 2 : coefficient of determination, RMSE: root mean square error, MAE: mean absolute error, PBIAS: percent bias, NSE: Nash–Sutcliffe efficiency coefficient, LCCC: Lin’s concordance correlation coefficient, and d: Willmott’s index of agreement. Based on F-test, at a probability of 0.05 (p< 0.05), significance of equation parameters for each response variable was determined. The high d(0.9620) values for ET a indicated that there was a good agreement between simulated and measured ET a . In general, the ET a model showed excellent agreement (LCCC 0.9), high performance, and low RMSE (0.498) (Table 3). The PBIAS and MAE values for ET a were between − 13.32% and 0.413, confirming the close agreement. Consequently, this result describes the ability to simplify and accurately predict ET a in deficit environments. Applications of regression model analysis with environmental variables are common in the literature in dry forests [ 118 – 120 ]. However, these implementations of variables in previous studies may be challenging to acquire for specific locations, requiring more simplified models. Our results also indicate the relationship between Hand the LST of ecosystems to determine ETa, an important variable in the energy balances of environments worldwide. Remote Sens. 2022,14, 1911 21 of 27 4. Conclusions In this study, we point out a significant tendency to increase the agricultural areas, which results in the progressive decrease of the Brazilian Caatinga biome. The vegetation cover is directly influenced by the soil–water regime; years of higher rainfall result in a lower percentage of suppression of the native forest in the municipality of Petrolina, Pernambuco (Brazil). The areas with pasture class presented hotspots due to degradative processes and higher surface temperatures, influenced by the sensible heat flux. A gradual increase in LST is observed in the municipality and it may cause future risks to forest areas. The SEBAL algorithm used in a semi-arid environment is a helpful tool to determine the energy and mass fluxes in different ecosystems. Notably, the Caatinga biome has particularities in biophysical parameters, according to the land cover and soil exposure on intra and inter-annual scales. The heterogeneity of the surface of the municipality of Petrolina, as a function of land use and land cover patterns, alters the energy exchange with the atmosphere. Our results also suggest a simplified and validated model for ET a determination in a semi-arid environment. The regression model could accurately predict the spatial distribution of ETa, with high R2and LCCC and low RMSE value. Thus, it is possible to suggest that the implementation of agricultural activities in the Petrolina should be carried out in a planned and sustainable way in order to mitigate the impacts that anthropic action causes on the Caatinga, especially with the increased vulnerability of this biome to the desertification process. However, further research is needed to investigate the spatial variations of the types of crops covering the soil in the municipality, as well as the dynamics of fires and their impacts on the diversity of the Caatinga biome. Field surveys and the use of unmanned aerial systems (UAS) could provide more detailed information at an intermediate and fine scale. Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/rs14081911/s1, Figure S1: Land use/land cover changes in Petrolina between 2013 and 2019. Positive values indicate an expansion of the respective land cover, negative values a contraction. The vertical axe is in million hectares (Mha). Please see Supplementary Table S1 to access values of land use/land cover changes; Table S1: Areas of expansion and contraction of land use/land cover changes in Petrolina between 2013 and 2019. Positive values indicate an expansion of the respective land cover, negative values a contraction. The values are in million hectares (Mha). Author Contributions: Conceptualization, A.M.d.R.F.J. and G.d.N.A.J.; methodology, M.V.d.S., A.d.S., and A.M.d.R.F.J.; software, M.V.d.S., A.d.S., and A.M.d.R.F.J.; validation, J.F.d.O.-J. and A.H.d.C.T.; investigation, J.L.B.d.S., H.P., P.E.T., L.S.B.d.S., and C.A.d.S.J.; data curation, M.V.d.S., A.d.S., A.M.d.R.F.J., and G.d.N.A.J.; writing—original draft preparation, A.M.d.R.F.J.; writing— review and editing, T.G.F.d.S., J.L.M.P.d.L., and E.A.S.; visualization, J.L.M.P.d.L. and A.H.d.C.T.; supervision, A.H.d.C.T., T.G.F.d.S., and J.F.d.O.-J.; project administration, T.G.F.d.S. and J.L.M.P.d.L.; funding acquisition, J.L.M.P.d.L. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by the Portuguese Foundation for Science and Technology (FCT), through projects ASHMOB (CENTRO-01-0145-FEDER-029351), GOLis (PDR2020-101-030913, Partnership nr. 344/Initiative nr. 21), MUSSELFLOW (PTDC/BIA-EVL/29199/2017), MEDWATERICE (PRIMA/0006/2018), and through the strategic project UIDB/04292/2020 granted to MARE—Marine and Environmental Sciences Centre, University of Coimbra, Coimbra, Portugal. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Landsat-8 image courtesy of the USGS/NASA (https://www.usgs. gov/core-science-systems/nli/landsat, accessed on 10 August 2021); MapBiomas data presented in this study are available at websites of Brazilian Annual Land Use and Land Cover Mapping Project (https://mapbiomas.org/en/project, accessed on 20 August 2021); and the meteorological data presented in this study are available at the website of the National Institute of Meteorology (https://portal.inmet.gov.br/, accessed on 10 August 2021). Remote Sens. 2022,14, 1911 22 of 27 Acknowledgments: The authors would like to thank the Research Support Foundation of the Pernambuco State (FACEPE, Brazil—APQ-0215-5.01/10 and FACEPE - APQ-1159-1.07/14), the National Council for Scientific and Technological Development (CNPq, Brazil) and also funds through the fellowship of the Research Productivity Program (CNPq 305286/2015-3, 304060/2016-0, 309681/20197, and 303767/2020-0), and the Coordination for the Improvement of Higher Education Personnel (CAPES, Brazil - Finance Code 001) for the research and study grants. The authors are also grateful for financial support from the Portuguese Foundation for Science and Technology (FCT), and the University of Coimbra, Portugal. In addition, we also would like to thank the anonymous reviewers for their insightful comments, of which significantly increased the value of this study. Conflicts of Interest: The authors declare no conflict of interest. Abbreviations Summary of all the symbols and acronyms used in this paper. Item Description aand bAre the calibration coefficients CpSpecific heat of air dWillmott’s index of agreement DEM Digital elevation model dT Near-surface air temperature gradient eaActual atmospheric water vapor pressure esWater vapor saturation pressure ETaActual evapotranspiration GSoil heat flux GEE Google Earth Engine HSensible heat flux HR Instantaneous relative humidity kvon Karman constant LAI Leaf area index LCCC Lin’s concordance correlation coefficient LE Latent heat flux LSE Land surface emissivity LST Land surface temperature LULC Land use and land cover MAE Mean absolute error NDVI Normalized Difference Vegetation Index NSE Nash-Sutcliffe efficiency coefficient PBIAS Percent bias PCA Principal component analysis rah Near-surface aerodynamic resistance to heat transport RMSE Root mean square error RnNet radiation Rn24 Daily net radiation R2Coefficient of determination SAVI Soil-Adjusted Vegetation Index T0Instantaneous air temperature TaAir temperature u* Friction velocity u200 Wind speed at the height of 200 m VCVegetation cover WPrecipitable water Remote Sens. 2022,14, 1911 23 of 27 Item Description z1and z2Are the two heights between the surface of the anchor pixels zom Momentum roughness length αsup Surface albedo εaAtmospheric emissivity ΛEvaporative fraction λLatent heat of vaporization of water ρair Air density ψmand ψhStability correction factors for momentum and sensible heat, respectively References 1. Arnan, X.; Leal, I.R.; Tabarelli, M.; Andrade, J.F.; Barros, M.F.; Câmara, T.; Jamelli, D.; Knoechelmann, C.M.; Menezes, T.G.C.; Menezes, A.G.S.; et al. A framework for deriving measures of chronic anthropogenic disturbance: Surrogate, direct, single and multi-metric indices in Brazilian Caatinga. Ecol. Indic. 2018,94, 274–282. [CrossRef] 2. Ferreira, T.R.; Silva, B.B.D.; De Moura, M.S.B.; Verhoef, A.; Nóbrega, R.L.B. The use of remote sensing for reliable estimation of net radiation and its components: A case study for contrasting land covers in an agricultural hotspot of the Brazilian semiarid region. Agric. For. Meteorol. 2020,291, 108052. [CrossRef] 3. Moro, M.F.; Nic Lughadha, E.; de Araújo, F.S.; Martins, F.R. A Phytogeographical Metaanalysis of the Semiarid Caatinga Domain in Brazil. Bot. Rev. 2016,82, 91–148. [CrossRef] 4. Magalhães, K.D.N.; Guarniz, W.A.S.; Sá, K.M.; Freire, A.B.; Monteiro, M.P.; Nojosa, R.T.; Bieski, I.G.C.; Custódio, J.B.; Balogun, S.O.; Bandeira, M.A.M. Medicinal plants of the Caatinga, northeastern Brazil: Ethnopharmacopeia (1980–1990) of the late professor Francisco Joséde Abreu Matos. J. Ethnopharmacol. 2019,237, 314–353. [CrossRef] 5. de Medeiros e Silva, É.; Paixão, V.H.F.; Torquato, J.L.; Lunardi, D.G.; de Oliveira Lunardi, V. Fruiting phenology and consumption of zoochoric fruits by wild vertebrates in a seasonally dry tropical forest in the Brazilian Caatinga. Acta Oecologica 2020 ,105, 103553. [CrossRef] 6. dos Santos, L.R.; Nascimento Lima, A.M.; Cunha, J.C.; Rodrigues, M.S.; Barros Soares, E.M.; dos Santos, L.P.A.; da Silva, A.V.L.; Ferreira Fontes, M.P. Does irrigated mango cultivation alter organic carbon stocks under fragile soils in semiarid climate? Sci. Hortic. 2019,255, 121–127. [CrossRef] 7. de Souza Leão, P.C.; do Nascimento, J.H.B.; de Moraes, D.S.; de Souza, E.R. Yield components of the new seedless table grape ‘BRS Ísis’ as affected by the rootstock under semi-arid tropical conditions. Sci. Hortic. 2020,263, 109114. [CrossRef] 8. Rallo, G.; Paço, T.A.; Paredes, P.; Puig-Sirera, À.; Massai, R.; Provenzano, G.; Pereira, L.S. Updated single and dual crop coefficients for tree and vine fruit crops. Agric. Water Manag. 2021,250, 106645. [CrossRef] 9. Gomes, L.S.; Maia, A.G.; de Medeiros, J.D.F. Fuzzified hedging rules for a reservoir in the Brazilian semiarid region. Environ. Chall. 2021,4, 100125. [CrossRef] 10. Marengo, J.A.; Torres, R.R.; Alves, L.M. Drought in Northeast Brazil—Past, present, and future. Theor. Appl. Climatol. 2017 ,129, 1189–1200. [CrossRef] 11. Martins, M.A.; Tomasella, J.; Rodriguez, D.A.; Alvalá, R.C.S.; Giarolla, A.; Garofolo, L.L.; Júnior, J.L.S.; Paolicchi, L.T.L.C.; Pinto, G.L.N. Improving drought management in the Brazilian semiarid through crop forecasting. Agric. Syst. 2018 ,160, 21–30. [CrossRef] 12. de Queiroz, M.G.; da Silva, T.G.F.; de Souza, C.A.A.; da Rosa Ferraz Jardim, A.M.; Araújo Júnior, G.D.N.; Souza, L.S.B.; Moura, M.S.B. Composition of Caatinga Species under Anthropic Disturbance and Its Correlation with Rainfall Partitioning. Floresta Ambient. 2021,28, 20190044. [CrossRef] 13. Barlow, J.; Lennox, G.D.; Ferreira, J.; Berenguer, E.; Lees, A.C.; Nally, R.M.; Thomson, J.R.; de Barros Ferraz, S.F.; Louzada, J.; Oliveira, V.H.F.; et al. Anthropogenic disturbance in tropical forests can double biodiversity loss from deforestation. Nature 2016 , 535, 144–147. [CrossRef] 14. da Silva Junior, C.A.; Teodoro, P.E.; Delgado, R.C.; Teodoro, L.P.R.; Lima, M.; de Andréa Pantaleão, A.; Baio, F.H.R.; De Azevedo, G.B.; de Oliveira Sousa Azevedo, G.T.; Capristo-Silva, G.F.; et al. Persistent fire foci in all biomes undermine the Paris Agreement in Brazil. Sci. Rep. 2020,10, 16246. [CrossRef] 15. da Rosa Ferraz Jardim, A.M.; da Silva, T.G.F.; de Souza, L.S.B.; do Nascimento Araújo Júnior, G.; Alves, H.K.M.N.; de SáSouza, M.; de Araújo, G.G.L.; de Moura, M.S.B. Intercropping forage cactus and sorghum in a semi-arid environment improves biological efficiency and competitive ability through interspecific complementarity. J. Arid Environ. 2021,188, 104464. [CrossRef] 16. Costa, M.D.S.; De Oliveira-Júnior, J.F.; Dos Santos, P.J.; Correia Filho, W.L.F.; De Gois, G.; Blanco, C.J.C.; Teodoro, P.E.; da Silva, C.A., Jr.; Santiago, D.D.B.; Souza, E.D.O.; et al. Rainfall extremes and drought in Northeast Brazil and its relationship with El Niño–Southern Oscillation. Int. J. Climatol. 2021,41, E2111–E2135. [CrossRef] 17. Tomasella, J.; Vieira, R.M.; Barbosa, A.A.; Rodriguez, D.A.; Santana, M.D.O.; Sestini, M.F. Desertification trends in the Northeast of Brazil over the period 2000–2016. Int. J. Appl. Earth Obs. Geoinf. 2018,73, 197–206. [CrossRef] 18. Vieira, R.M.S.P.; Tomasella, J.; Alvalá, R.C.S.; Sestini, M.F.; Affonso, A.G.; Rodriguez, D.A.; Barbosa, A.A.; Cunha, A.P.M.A.; Valles, G.F.; Crepani, E.; et al. Identifying areas susceptible to desertification in the Brazilian northeast. Solid Earth 2015 ,6, 347–360. [CrossRef] Remote Sens. 2022,14, 1911 24 of 27 19. Ribeiro, K.; de Sousa-Neto, E.R.; de Carvalho, J.A.; Lima, J.R.D.S.; Menezes, R.; Duarte-Neto, P.J.; Guerra, G.D.S.; Ometto, J.P.H.B. Land cover changes and greenhouse gas emissions in two different soil covers in the Brazilian Caatinga. Sci. Total Environ. 2016 , 571, 1048–1057. [CrossRef] 20. de Castro Teixeira, A.H.; Leivas, J.F.; Andrade, R.G.; Hernandez, F.B.T. Water productivity assessments with Landsat 8 images in the Nilo Coelho irrigation scheme. IRRIGA 2015,1, 1–10. [CrossRef] 21. Ronquim, C.C.; Leivas, J.F.; de Castro Teixeira, A.H.; Silva, G.B.; Garçon, E.A.M. Water indicators based on SPOT 6 satellite images in irrigated area at the Paracatu River Basin, Brazil. In Remote Sensing for Agriculture, Ecosystems, and Hydrology XIX, 104211I; International Society for Optics and Photonics: Bellingham, WA, USA, 2017; Volume 10421, p. 104211I. 22. Cunha, J.; Nóbrega, R.L.B.; Rufino, I.; Erasmi, S.; Galvão, C.; Valente, F. Surface albedo as a proxy for land-cover clearing in seasonally dry forests: Evidence from the Brazilian Caatinga. Remote Sens. Environ. 2020,238, 111250. [CrossRef] 23. Barbosa, H.A.; Lakshmi Kumar, T.V.; Paredes, F.; Elliott, S.; Ayuga, J.G. Assessment of Caatinga response to drought using Meteosat-SEVIRI Normalized Difference Vegetation Index (2008–2016). ISPRS J. Photogramm. Remote Sens. 2019 ,148, 235–252. [CrossRef] 24. de Queiroz, M.G.; da Silva, T.G.F.; Zolnier, S.; da Rosa Ferraz Jardim, A.M.; de Souza, C.A.A.; do Nascimento Araújo Júnior, G.; de Morais, J.E.F.; de Souza, L.S.B. Spatial and temporal dynamics of soil moisture for surfaces with a change in land use in the semi-arid region of Brazil. Catena 2020,188, 104457. [CrossRef] 25. Fisher, J.B.; Melton, F.; Middleton, E.; Hain, C.; Anderson, M.; Allen, R.; McCabe, M.F.; Hook, S.; Baldocchi, D.; Townsend, P.A.; et al. The future of evapotranspiration: Global requirements for ecosystem functioning, carbon and climate feedbacks, agricultural management, and water resources. Water Resour. Res. 2017,53, 2618–2626. [CrossRef] 26. Souza, R.; Hartzell, S.; Feng, X.; Antonino, A.C.D.; de Souza, E.S.; Menezes, R.S.C.; Porporato, A. Optimal management of cattle grazing in a seasonally dry tropical forest ecosystem under rainfall fluctuations. J. Hydrol. 2020,588, 125102. [CrossRef] 27. Fendrich, A.N.; Barretto, A.; de Faria, V.G.; de Bastiani, F.; Tenneson, K.; Guedes Pinto, L.F.; Sparovek, G. Disclosing contrasting scenarios for future land cover in Brazil: Results from a high-resolution spatiotemporal model. Sci. Total Environ. 2020 ,742, 140477. [CrossRef] 28. Lopes, V.C.; Parente, L.L.; Baumann, L.R.F.; Miziara, F.; Ferreira, L.G. Land-use dynamics in a Brazilian agricultural frontier region, 1985–2017. Land Use Policy 2020,97, 104740. [CrossRef] 29. Blondeel, H.; Landuyt, D.; Vangansbeke, P.; De Frenne, P.; Verheyen, K.; Perring, M.P. The need for an understory decision support system for temperate deciduous forest management. For. Ecol. Manag. 2021,480, 118634. [CrossRef] 30. de Araujo, H.F.; Machado, C.C.; Pareyn, F.G.; Nascimento, N.F.D.; Araújo, L.D.; de A. P. Borges, L.A.; Santos, B.A.; Beirigo, R.M.; Vasconcellos, A.; Dias, B.D.O.; et al. A sustainable agricultural landscape model for tropical drylands. Land Use Policy 2021 ,100, 104913. [CrossRef] 31. Liu, S.; Su, H.; Zhang, R.; Tian, J.; Chen, S.; Wang, W. Regional Estimation of Remotely Sensed Evapotranspiration Using the Surface Energy Balance-Advection (SEB-A) Method. Remote Sens. 2016,8, 644. [CrossRef] 32. Mutti, P.R.; da Silva, L.L.; Medeiros, S.D.S.; Dubreuil, V.; Mendes, K.R.; Marques, T.V.; Lúcio, P.S.; e Silva, C.M.S.; Bezerra, B.G. Basin scale rainfall-evapotranspiration dynamics in a tropical semiarid environment during dry and wet years. Int. J. Appl. Earth Obs. Geoinf. 2019,75, 29–43. [CrossRef] 33. Teixeira, A.D.C.; de Miranda, F.; Leivas, J.; Pacheco, E.; Garçon, E. Water productivity assessments for dwarf coconut by using Landsat 8 images and agrometeorological data. ISPRS J. Photogramm. Remote Sens. 2019,155, 150–158. [CrossRef] 34. Moreira, E.B.M.; Nóbrega, R.S.; Da Silva, B.B.; Ribeiro, E.P. Estimativa da evapotranspiração em área urbana através de imagens digitais TM-Landsat 5. Geosul 2019,34, 559–585. [CrossRef] 35. da Silva, M.V.; Pandorfi, H.; de Almeida, G.L.P.; de Lima, R.P.; dos Santos, A.; da Rosa Ferraz Jardim, A.M.; Rolim, M.M.; da Silva, J.L.B.; Batista, P.H.D.; da Silva, R.A.B.; et al. Spatio-temporal monitoring of soil and plant indicators under forage cactus cultivation by geoprocessing in Brazilian semi-arid region. J. S. Am. Earth Sci. 2021,107, 103155. [CrossRef] 36. Mhawej, M.; Faour, G. Open-source Google Earth Engine 30-m evapotranspiration rates retrieval: The SEBALIGEE system. Environ. Model. Softw. 2020,133, 104845. [CrossRef] 37. Júnior, J.B.C.; e Silva, C.M.S.; De Almeida, H.A.; Bezerra, B.; Spyrides, M.H.C. Detecting linear trend of reference evapotranspiration in irrigated farming areas in Brazil’s semiarid region. Theor. Appl. Climatol. 2019,138, 215–225. [CrossRef] 38. IBGE Instituto Brasileiro de Geografia e Estatística. Available online: https://cidades.ibge.gov.br/brasil/pe/petrolina/panorama (accessed on 29 August 2021). 39. NASA Giovanni. National Aeronautics and Space Administration. Available online: https://giovanni.gsfc.nasa.gov/giovanni/ (accessed on 29 August 2021). 40. Santos, C.; da Silva, R.M.; Silva, A.M.; Neto, R.M.B. Estimation of evapotranspiration for different land covers in a Brazilian semi-arid region: A case study of the Brígida River basin, Brazil. J. S. Am. Earth Sci. 2017,74, 54–66. [CrossRef] 41. Tasumi, M. Progress in Operational Estimation of Regional Evapotranspiration Using Satellite Imagery; University of Idaho: Moscow, ID, USA, 2003. 42. Moletto-Lobos, I.; Mattar, C.; Barichivich, J. Performance of Satellite-Based Evapotranspiration Models in Temperate Pastures of Southern Chile. Water 2020,12, 3587. [CrossRef] 43. Consoli, S.; Inglese, P.; Inglese, G. Determination of evapotranspiration and crop coefficient of cactus pear (Opuntia ficus-indica Mill.) with an energy balance technique. Acta Hortic. 2013,995, 117–124. [CrossRef] Remote Sens. 2022,14, 1911 25 of 27 44. Liu, J.; You, Y.; Li, J.; Sitch, S.; Gu, X.; Nabel, J.E.M.S.; Lombardozzi, D.; Luo, M.; Feng, X.; Arneth, A.; et al. Response of global land evapotranspiration to climate change, elevated CO 2 , and land use change. Agric. For. Meteorol. 2021 ,311, 108663. [CrossRef] 45. Hartzell, S.; Bartlett, M.S.; Porporato, A. Unified representation of the C3, C4, and CAM photosynthetic pathways with the Photo3 model. Ecol. Model. 2018,384, 173–187. [CrossRef] 46. Laipelt, L.; Ruhoff, A.L.; Fleischmann, A.; Kayser, R.H.B.; Kich, E.D.M.; Da Rocha, H.R.; Neale, C.M.U. Assessment of an Automated Calibration of the SEBAL Algorithm to Estimate Dry-Season Surface-Energy Partitioning in a Forest–Savanna Transition in Brazil. Remote Sens. 2020,12, 1108. [CrossRef] 47. Allen, R.; Waters, R.; Bastiaanssen, W.; Tasumi, M.; Trezza, R. SEBAL (Surface Energy Balance Algorithms for Land)—Idaho Implementation, Advanced Training and Users Manual, Version 1.0; Idaho Department of Water Resources: Boise, ID, USA, 2002. 48. Bastiaanssen, W.G.M. SEBAL-based sensible and latent heat fluxes in the irrigated Gediz Basin, Turkey. J. Hydrol. 2000 ,229, 87–100. [CrossRef] 49. Bastiaanssen, W.G.M.; Noordman, E.J.M.; Pelgrum, H.; Davids, G.; Thoreson, B.P.; Allen, R.G. SEBAL Model with Remotely Sensed Data to Improve Water-Resources Management under Actual Field Conditions. J. Irrig. Drain. Eng. 2005 ,131, 85–93. [CrossRef] 50. Allen, R.G.; Tasumi, M.; Trezza, R. Satellite-Based Energy Balance for Mapping Evapotranspiration with Internalized Calibration (METRIC)—Model. J. Irrig. Drain. Eng. 2007,133, 380–394. [CrossRef] 51. Cheng, M.; Jiao, X.; Li, B.; Yu, X.; Shao, M.; Jin, X. Long time series of daily evapotranspiration in China based on the SEBAL model and multisource images and validation. Earth Syst. Sci. Data 2021,13, 3995–4017. [CrossRef] 52. Filho, W.L.F.C.; Santiago, D.D.B.; de Oliveira-Júnior, J.F.; Junior, C.A.D.S. Impact of urban decadal advance on land use and land cover and surface temperature in the city of Maceió, Brazil. Land Use Policy 2019,87, 104026. [CrossRef] 53. Rouse, J.W.; Haas, R.H.; Schell, J.A.; Deering, D.W.; Harlan, J.C. Monitoring the Vernal Advancement and Retrogradation (Greenwave Effect) of Natural Vegetation. NASA/GSFCT Type III Final Report; NASA/GSFCT: Greenbelt, MD, USA, 1974; pp. 1–390. 54. Gao, Q.; Li, Y.; Wan, Y.; Lin, E.; Xiong, W.; Jiangcun, W.; Wang, B.; Li, W. Grassland degradation in Northern Tibet based on remote sensing data. J. Geogr. Sci. 2006,16, 165–173. [CrossRef] 55. de Lima, I.P.; Jorge, R.G.; de Lima, J.L.M.P. Remote Sensing Monitoring of Rice Fields: Towards Assessing Water Saving Irrigation Management Practices. Front. Remote Sens. 2021,2, 762093. [CrossRef] 56. Teixeira, A.H.C.; Bastiaanssen, W.G.M.; Ahmad, M.D.; Bos, M.G. Reviewing SEBAL input parameters for assessing evapotranspiration and water productivity for the Low-Middle São Francisco River basin, Brazil: Part B: Application to the regional scale. Agric. For. Meteorol. 2009,149, 477–490. [CrossRef] 57. Bright, R.M.; Davin, E.; O’Halloran, T.; Pongratz, J.; Zhao, K.; Cescatti, A. Local temperature response to land cover and management change driven by non-radiative processes. Nat. Clim. Chang. 2017,7, 296–302. [CrossRef] 58. Bastiaanssen, W.G.M.; Pelgrum, H.; Soppe, R.W.O.; Thoreson, B.P.; Allen, R.G.; Teixeira, A.H.C. Thermal-infrared technology for local and regional scale irrigation analyses in horticultural systems. Acta Hortic. 2008,792, 33–46. [CrossRef] 59. Teixeira, A.H.C.; Bastiaanssen, W.G.M.; Moura, M.S.B.; Soares, J.M.; Ahmad, M.D.; Bos, M.G. Energy and water balance measurements for water productivity analysis in irrigated mango trees, Northeast Brazil. Agric. For. Meteorol. 2008 ,148, 1524–1537. [CrossRef] 60. Filho, W.L.F.C.; De Oliveira-Júnior, J.F.; De Barros Santiago, D.; De Bodas Terassi, P.M.; Teodoro, P.E.; De Gois, G.; Blanco, C.J.C.; De Almeida Souza, P.H.; da Silva Costa, M.; Gomes, H.B.; et al. Rainfall variability in the Brazilian northeast biomes and their interactions with meteorological systems and ENSO via CHELSA product. Big Earth Data 2019,3, 315–337. [CrossRef] 61. Teixeira, A.D.C.; Bastiaanssen, W.; Ahmad, M.-U.; Bos, M. Reviewing SEBAL input parameters for assessing evapotranspiration and water productivity for the Low-Middle São Francisco River basin, Brazil: Part A: Calibration and validation. Agric. For. Meteorol. 2009,149, 462–476. [CrossRef] 62. Alvares, C.A.; Stape, J.L.; Sentelhas, P.C.; de Moraes Gonçalves, J.L.; Sparovek, G. Köppen’s climate classification map for Brazil. Meteorol. Z. 2013,22, 711–728. [CrossRef] 63. Beck, H.E.; Zimmermann, N.E.; McVicar, T.R.; Vergopolan, N.; Berg, A.; Wood, E.F. Present and future Köppen-Geiger climate classification maps at 1-km resolution. Sci. Data 2018,5, 180214. [CrossRef] 64. Oliveira, P.T.; e Silva, C.M.S.; Lima, K.C. Climatology and trend analysis of extreme precipitation in subregions of Northeast Brazil. Theor. Appl. Climatol. 2017,130, 77–90. [CrossRef] 65. da Rosa Ferraz Jardim, A.M.; da Silva, M.V.; Silva, A.R.; dos Santos, A.; Pandorfi, H.; de Oliveira-Júnior, J.F.; de Lima, J.L.; de Souza, L.S.B.; do Nascimento Araújo Júnior, G.; Lopes, P.M.O.; et al. Spatiotemporal climatic analysis in Pernambuco State, Northeast Brazil. J. Atmos. Solar-Terr. Phys. 2021,223, 105733. [CrossRef] 66. Preston, W.; Nascimento, C.; Silva, Y.; Silva, D.J.; Ferreira, H.A. Soil fertility changes in vineyards of a semiarid region in Brazil. J. Soil Sci. Plant Nutr. 2017,17, 672–685. [CrossRef] 67. Menezes, K.M.S.; Silva, D.K.A.; Gouveia, G.V.; da Costa, M.M.; Queiroz, M.A.A.; Yano-Melo, A.M. Shading and intercropping with buffelgrass pasture affect soil biological properties in the Brazilian semi-arid region. Catena 2019,175, 236–250. [CrossRef] 68. Giongo, V.; Coleman, K.; da Silva Santana, M.; Salviano, A.M.; Olszveski, N.; Silva, D.J.; Cunha, T.J.F.; Parente, A.; Whitmore, A.P.; Richter, G.M. Optimizing multifunctional agroecosystems in irrigated dryland agriculture to restore soil carbon—Experiments and modelling. Sci. Total Environ. 2020,725, 138072. [CrossRef] 69. INMET Instituto Nacional de Meteorologia. Available online: https://portal.inmet.gov.br/ (accessed on 29 August 2021).