Climatic impacts of bushland to cropland conversion in Eastern Africa

• Climatic impacts of bushland loss estimated using satellite image and


Introduction
Land cover change has long been recognized as a key driver of anthropogenic climate change (Ciais et al., 2013;Bonan, 2008). According to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC), land cover change was identified as the second biggest driver for increased carbon dioxide (CO 2 ) emissions next to the burning of fossil fuels (Ciais et al., 2013). Land cover change affects climate through the modification of surface biophysical properties that have a direct influence on the exchange of water, energy, and carbon between the land surface and the atmosphere (Ciais et al., 2013;Davin and de Noblet-Ducoudre, 2010;Eltahir and Bras, 1966). A recent observational study reported a global average increase in land surface temperature of 0.23 ± 0.03°C following land cover changes during the period 2000 to 2015; the warming was mainly attributed to agricultural expansion in the tropics (Duveiller et al., 2018).
Tropical regions have undergone considerable changes in woody vegetation cover due to agricultural expansion in recent decades (Brink et al., 2014;Pfeifer et al., 2012). For instance, between 2002 and2008 in East Africa, woody vegetation cover decreased by 5. 1%, 15.8%, and 19.4% from forests, woodland, and shrubland, respectively (Pfeifer et al., 2012). This change is even more pronounced in the Horn of Africa, where croplands increased by 28% at the expense of natural woody vegetation between 1990 and 2010; the annual rate was 1.4% (Brink et al., 2014). The Horn of Africa is mainly covered by deciduous bushland and thicket ecoregion that consists of Acacia and Commiphora trees (Bouvet et al., 2018;Dinerstein et al., 2017;White, 1983) (Fig. 1a). This ecoregion constitutes much of the woody cover, since other woody vegetation classes (such as forest and savannas) cover smaller areas in the highlands and mountains. The woody vegetation cover in this region is declining particularly due to bushland clearing associated with strong anthropogenic pressure from cropland encroachment (Pellikka et al., 2013).
However, the impacts of bushland-to-cropland conversion on surface biophysical properties and their consequent effect on surface energy balance and land surface temperature (LST) is currently Fig. 1. Spatial extents of (a) Acacia-Commiphora bushlands and thickets and Montane forest in East Africa, modified from Dinerstein et al. (2017) and based on White (1983). Panel (b) shows the location of the study areas (Mwatate and Ndome). The geographic extent of panel (b) shows the area of interest (AOI) used for evapotranspiration (ET) modeling, from which ET values for the study areas were extracted. A 30-m digital surface model (DSM) from the Japan Aerospace Exploration Agency (JAXA) was used as background. Panels (c) and (d) show bushland and cropland distribution in Mwatate and Ndome as derived from RapidEye imagery. Masked pixels represent unstable areas and slope N 5°. Panel (e) displays bushlands cleared by small-scale farmers for subsistence agriculture in Mwatate with Chawia Mountain in background, and panel (f) shows Acacia-Commiphora bushland in Mwatate. Photos: Professor Petri Pellikka, 2019. uncertain. Many reasons contribute to this knowledge gap. First, previous (observational and modeling) studies at global and regional scale were too coarse (≥1 km; Abera et al., 2019;Duveiller et al., 2018;Bright et al., 2017) to study the specific role of bushlands. Second, croplands often occur mixed with other vegetation classes in the region at this resolution (Zhang et al., 2005), making the computation of biophysical effects of land cover change inconsistent except for the conversion of forest to cropland and savanna to grassland (Abera et al., 2019). Third, the lack of high (spatial and temporal) resolution biophysical and energy flux data further contributes to this knowledge gap.
Recent advances in satellite observation and modeling present new opportunities to quantify biophysical (e.g., albedo and evapotranspiration) and temperature changes associated with land cover transitions in data-scarce regions at more detailed scales than before. For instance, studies in the integration of Landsat surface reflectance and MODIS BRDF anisotropy information able to generate land surface albedo at fine spatial resolution (30-m) (Shuai et al., 2011). In addition, progress in remote-sensing-based energy balance models enable estimation of changes in actual evapotranspiration at higher resolution over large areas (Allen et al., 2007). Mapping Evapotranspiration at High Resolution with Internalized Calibration (METRIC) is one of the remotesensing based, single-source energy balance models that have produced satisfactory accuracy in different continents and complicated terrains (refer McShane et al., 2017 for detailed review and evaluation of this model). Furthermore, the capacity of METRIC to solve all energy balance terms in tandem with its limited requirement for the number of weather stations (i.e., only 1 to 2 weather stations per Landsat scene; Allen et al., 2007) makes it useful in data-scarce regions, such as the Horn of Africa. Therefore, through harnessing the recent advances in remote sensing and modeling, further studies at smaller scale are possible in areas that have undergone considerable conversion from bushland to cropland to identify their biophysical effects on local climate.
In this study, we used high-resolution satellite imagery and land surface flux modeling with the objective of identifying the biophysical impacts of bushland clearing on surface energy balance and land surface temperature. We proposed the hypothesis that bushland-to-cropland conversion affects local surface energy balance and induces significant warming. The hypothesis was tested in southern Kenya through conducting detailed analysis on (1) how the actual (blue-sky) shortwave albedo changes during bushland-to-cropland conversion using Landsat 8 and MODIS BRDF anisotropy data, (2) how the actual evapotranspiration changes during this conversion, and the associated surface energy flux and LST change.

Study sites
The study sites (Mwatate and Ndome) are located in the lowlands of Taita Hills, southern Kenya, and cover an area of 50 km 2 and 37 km 2 , respectively ( Fig. 1). They were selected based on the prominent occurrence of croplands and bushlands, which are not mixed with other land use classes (e.g., settlement). The lowlands are characterized by a semi-arid climate and receive a mean annual precipitation between 300 and 600 mm with an average temperature of 22°C (Bouvet et al., 2018). The dry seasons are hottest during January and February and cooler from June to October ( Supplementary Fig. S1a). Vegetation in the lowlands is characterized by bushlands, thickets, and croplands . The bushlands and thickets belong to the Somali-Masai Acacia-Commiphora deciduous bushland and thicket ecoregion White, 1983). The bushlands, which are the main vegetation type in the study areas, form an open stand with a dense canopy cover of 40% or more and are 3 to 5 m tall (White, 1983;Pellikka et al., 2018; Fig. 1e, f).

Remote sensing and weather station data
We used satellite observation data to derive biophysical changes (i.e., albedo and ET changes) and their associated impact on climate (i.e., LST) due to bushland-to-cropland conversion. Landsat 8, MODIS (BRDF and Aerosol Optical Depth), RapidEye, and JAXA digital surface model (DSM) were used in this study.
Landsat 8 Operational Land Imager (OLI) provides multispectral imagery at 30-m spatial and 16-day temporal resolutions. Thermal Infrared Sensor (TIRS) provides additional thermal bands at 100-m resolution (distributed at 30-m resolution corresponding to OLI bands). In this study, Landsat 8 Collection 1 level 1, level 2 (surface reflectance), and thermal band were used to compute ET in METRIC modeling. The surface reflectance data were also used in the blue-sky albedo calculation in combination with MODIS BRDF. Furthermore, band 10 of TIRS was employed in LST retrieval. For these purposes in total, eight nearly cloud-free (b10%) Landsat 8 scenes were downloaded (https://earthexplorer.usgs.gov/) and used during 2014 to 2018. While a detailed analysis was performed by applying the first scene (2014-11-27), the remaining was used to verify the consistency of the result from multiple dates.
MODIS MCD43A1 daily product, at 500-m resolution, provides the model weighting parameters (isotropic, volumetric, and geometric) required to derive the BRDF and albedo using surface reflectance as input (https://lpdaac.usgs.gov/products/mcd43a1v006/). This product was used to calculate the blue-sky albedo together with Landsat 8 surface reflectance and Aerosol Optical Depth (AOD) data. The AOD data were obtained from the MCD19A2 daily product at 0.55 μm and 1 km (https:// lpdaac.usgs.gov/products/mcd19a2v006/). It is a Multi-Angle Implementation of Atmospheric Correction (MAIAC) algorithm-based Level-2 gridded AOD data, derived from combined Terra and Aqua MODIS (Lyapustin and Wang, 2018). The MODIS products (MCD43A1 and AOD) were resampled to 300-m to match our analysis grid (see Sections 2.3.1 and 2.3.4).
For identifying bushland and cropland areas, we used a highresolution (5 m) RapidEye multi-spectral image that provides bands in the blue (440-510 nm), green (520-590 nm), red (630-685 nm), red edge (690-730 nm), and near infrared (760-850 nm) wavelength. The RapidEye image was obtained from https://www.planet.com/ (Planet Team, 2017). Prior to use, an atmospheric correction was applied using the Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) algorithm. A digital surface model (DSM) from JAXA was used for parameterizing the METRIC model and for topographic correction of Landsat 8 surface reflectance data. This 30-m resolution DSM is a gap-filled product and is freely available from http://www.eorc.jaxa.jp/ALOS/en/aw3d30/.
Two weather stations data that fall in the same Landsat scene (Path 167 and row 62) were used; these were located in Maktau and in Voi town (Fig. 1b). The Maktau station, which is located in a maize field, and Voi station are operated by the Taita Research Station of the University of Helsinki and the Kenya Meteorological Department, respectively. The Maktau weather station data are available at 30-min intervals (2014-2018) and hence were used for METRIC modeling. While the daily Voi weather data was applied to calculate actual ET by FAO-PM (Penman Monteith) equations, its result was used for comparison with the METRIC daily ET.

Identifying adjacent bushland and cropland
First, to identify cropland and bushland across the study area, a supervised image segmentation and classification was performed using a high-resolution (5-m) RapidEye daily image based on prior knowledge of the study area. Accuracy of the classification was assessed using a systematic sample taken from Google Earth image (100 samples for each study sites; 200 samples in total). Overall thematic accuracy was 88% for Mwatate and 89% for Ndome.
The result was aggregated to 30-m to match Landsat 8 resolution. The fractions of cropland and bushland were then calculated for 300 × 300 m pixels across the study area (Fig. 2d,i). All the pixels with a fraction of cropland/bushland b25% were discarded to minimize impact of classification accuracy on results and assure the co-occurrence of the both classes in each pixel (Fig. 2e, j). The temporal stability of pixels was further checked on a pixel-by-pixel basis using a Google Earth image across 2014 to 2018 to ensure that only pixels without land cover change in the studied period are selected. By this process, a total of 11 pixels that have undergone land cover change were discarded from Mwatate and Ndome (black pixels in Fig. 2e, j). Finally, as topography can have an impact on the analysis, we removed all pixels having an average slope N 5° (Fig. 2e). Because of this, the mean elevation differences between cropland and shrubland within a 300 × 300 m grid were constrained to b20 m ( Supplementary  Fig. S1b). Furthermore, pixels affected by clouds and cloud shadows were masked using pixel quality assurance band prior to analysis from each daily Landsat 8 scenes between 2014 and 2018. In total, 28,010 pixels in Mwatate and 12,750 pixels in Ndome were masked from these scenes. For detecting cloud and cloud shadows, Landsat 8 uses an improved Fmask algorithm (LASRC, Landsat 8 surface reflectance code product guide, 2019). For details of the algorithm refer Zhu et al. (2015).

Estimating blue-sky albedo
The Landsat 8 blue-sky visible (VIS), near-infrared (NIR), and shortwave (SW) broadband albedo was calculated from the MODIS BRDF anisotropy and Landsat 8 surface reflectance data (Shuai et al., 2011). The absolute accuracy of applying this approach can range from 2% to 5% compared with field measurement (Shuai et al., 2011) The MODIS BRDF data was used to calculate the albedo to reflectance ratio (ATRR) at Landsat 8 view geometry (i.e., sensor and sun view angle). For this task, first the white-sky and black-sky albedos were calculated for each of the Landsat 8 equivalent MODIS bands (Supplementary  Table S1). Then the MODIS reflectance at Landsat view geometry was obtained using: where R is MODIS reflectance; ʎ is wavelength; θ v is sensor view angle; θ s is solar zenith angle; φ is the relative azimuth; f iso , f vol , and f geo are the isometric, volumetric, and geometric scattering parameters from MODIS BRDF, respectively; and K vol and K geo are the coefficient of volumetric and geometric kernels calculated at Landsat 8 view geometries, respectively (Li and Strahler, 1992;Ross, 1981). Second, the band-specific ATRR was applied to convert the Landsat 8 surface reflectance to black-sky and white-sky albedos. The blue-sky Landsat 8 albedo was computed from a linear combination of the black-sky (α bs ) and white-sky albedos (α ws ) as a function of diffuse skylight (Eq. (2)). The fraction of diffuse skylight (S) was obtained from the MODIS lookup table using MODIS Aerosol Optical Thickness at 550 nm and Landsat 8 solar zenith angle (θ).
Finally, the Landsat narrow bands were converted to broadband VIS, NIR, and SW blue-sky albedos (Liang, 2001) and used to estimate albedo changes during bushland-to-cropland conversion on 2014-11-27. See  -11-27 (b, g), and maps of croplands (yellow pixels) and bushland (green pixels) (c, h) in the study areas (Mwatate and Ndome) derived from RapidEye image, percentage of cropland within 300 × 300 m pixel across the study areas (d, i), and selected samples for analysis (e, j). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) spatial distribution of blue-sky albedo in the study areas in the Supplementary materials (Fig. S2a, b).

Estimating evapotranspiration using METRIC model
A field-scale ETa was estimated using a METRIC model (Allen et al., 2007) with an automated anchor pixel-selection procedure applied in Water package in R (Olmedo et al., 2016). METRIC is a single-source energy balance model that estimates ET as a residual of surface energy balance (Eqs. (3)- (6)). The energy balance is internally calibrated using a ground-based reference ET from a weather station within or near the study area for a Landsat scene (see Allen et al. (2007) for details of the algorithm).
where LE is the latent heat flux (W m −2 ); Rn is the surface net radiation flux (W m −2 ); G is the soil heat flux (W m −2 ); and H is the sensible heat flux (W m −2 ). Rn, G, and H were computed using weather station data, digital elevation model, and Landsat 8 as described in Allen et al. (2007). For example, G was calculated by multiplying G/Rn ratio (derived using empirical equation from Bastiaanssen (2000)) by Rn. The instantaneous evapotranspiration (ET inst ) at satellite overpass is computed for every pixel using: where ʎ is the latent heat of vaporization (ʎ = [2.501-0.00236 (Ts −273.15)] × 10 6 ; Ts is surface temperature); ρ w is the density of water (1000 kg m −3 ). The reference ET fraction (ET rf ) was computed from the ratio of ET inst and reference ET (ET r ) using Eq. (5). Reference ET was calculated from the Maktau weather station data (air-temperature, solar radiation, wind speed, relative humidity, location and elevation of the weather station, and height of wind speed measurement) applying standardized Penman-Monteith equation (ASCE-EWRI, 2005). The Maktau weather station is located in maize field.
As METRIC requires cloud-free scenes, all nearly cloud-free (b10% cloud cover) Landsat 8 scenes acquired between 2014 and 2018 were selected. Eight scenes were identified by this process. As all images need to be clear sky, all clouds, cloud shadows, and other features that might affect the ET modeling (cirrus, terrain occlusion, and water bodies) were then masked using the quality assurance band.
To avoid any translation of topographic effect (slope and aspect) to the two study areas, which are located in flat terrain, topographic correction was subsequently applied to the mountainous part of the METRIC AOI (Fig. 1b) using C-correction (Teillet et al., 1982) on Landsat 8 surface reflectance data.
LST was calculated from the first thermal band (band 10) of Landsat 8 TIRS using a single-channel algorithm (Jimenez-Munoz et al., 2009). Although the stray light problem that had affected the thermal detector was already corrected in Landsat 8 (i.e., the correction reduced the TIRS uncertainty b0.5%; Gerace and Montanaro, 2017), we preferred singlechannel (SC) algorithms as the use of band 11 for split-window algorithm was not yet recommended by the Landsat team until further testing (García-Santos et al., 2018). Emissivity value used in the computation of LST was derived as a function of LAI (Ɛ = 0.95 + 0.01 × LAI for LAI ≤ 3) as described in Tasumi (2003). Since LAI in the study areas are b3, we did not assume a constant emissivity (i.e. 0.98 for LAI N 3 as described in Tasumi, 2003). The LST, calculated with the SC algorithm, was also corrected for topographic effect using a standard lapse rate of 0.0065°C m −1 in METRIC modeling for the METRIC AOI (Fig. 1b) to reduce the impacts of temperature change associated with elevation difference.
The broadband albedo here (in METRIC) was estimated from the surface reflectance narrow bands using the Olmedo coefficient, which was specifically calculated for Landsat 8 using the Simple Model of the Atmospheric Radiative Transfer of Sunshine (SMARTS2) version 2.9.5 in Water Package (Olmedo et al., 2016) to maintain consistency in the use of the model.
Momentum Roughness Length (Z om = a1(NDVI/α) + b1), which is required to calculate the aerodynamic resistance (R ah ), was custom calculated from NDVI, albedo (α), surface model (slope and aspect), and regression constants (a1 and b1) derived from regression between ln (Z om ) and NDVI/α (Allen et al., 2007). The customization was required as our study contains bushlands, whose heights reach 3 to 5 m, and the original formulation in METRIC (Z om = 0.018LAI; Tasumi, 2003) is applicable for short agricultural crops (b1 m).
The leaf area index (LAI) was estimated from the METRIC 2010 method. This method estimates LAI from soil adjusted vegetation index (SAVI; Huete, 1988): where ρNIR and ρR are near infrared and red reflectance and L is soil correction factor, which was set to 0.5 for our study area as it is partly vegetated (Allen et al., 2007). Comparison of mean LAI estimate from METRIC (1.46 ± SD 0.74) against field measurements from previous study (1.28 ± SD 0.55; Itkonen, 2012) showed small difference over bushland (e.g. in Ndome). Furthermore, the impacts of LAI estimation error on METRIC energy balance calculation are small (Allen et al., 2007). In Water package, the cold and hot pixels are selected automatically based on defined criteria from Allen et al. (2011) andTasumi (2003) ( Table 1). Cold/hot pixels that do not meet these criteria will not be selected. This minimizes the uncertainty associated in the selection of cold and hot pixels as the process does not depend on experience of the user. The expected error in using automated method is around 10% compared to manual approach by experienced user in semi-arid environment . Furthermore, terrain slope and aspect from the digital elevation model were also used in the computation of other parameters that are used for ET modeling in METRIC (e.g., incoming shortwave, momentum roughness length, sensible heat flux, 24-hour ET) as described in Allen et al. (2007). See spatial distribution of ET in the study areas in the Supplementary materials (Fig. S3).
To check the robustness of the method, METRIC was compared with SEBAL model (Fig. S4). The comparison showed moderate agreement (r 2 = 0.65) and hence robustness of the method.

Assessing biophysical changes and impact on land surface temperature
The biophysical changes (albedo and ET changes) and the associated impact on LST was assessed for each pixel by calculating the mean difference in albedo, ET, and LST between cropland and bushland in each 300 × 300 m pixel using a space-for-time approach. The statistical significance of mean differences was tested using Welch's t-test.
Space-for-time approach identifies the potential impacts of land cover conversion between contrasting vegetation types that are located close to each other and share similar environmental conditions (Duveiller et al., 2018). The approach has been successfully applied at the global (e.g., Duveiller et al., 2018;Li et al., 2016;Li et al., 2015) and local level (e.g., Peng et al., 2014;Zhao and Jackson, 2017). Furthermore, it is more commonly used than the observed land cover change approach, which is limited by the accuracy of land cover map, frequency and extent of changes, and interannual climate variability signal (Duveiller et al., 2018). Its applicability for environmental studies has been also tested (Blois et al., 2013). However, there are some drawback to this approach. For example, it shows potential impact instead of actual impact.
As albedo and ET change affect the energy balance, the mean difference in the energy fluxes (shortwave, longwave, sensible, latent energy, and ground heat fluxes) were computed for the study areas. The relationship of LST changes with ET and albedo changes were also explored using a partial correlation and multiple regression approach. Furthermore, as ET and albedo change cannot fully explain the LST changes, the relative contribution of each of the energy balance terms on the LST changes were assessed using an energy balance decomposition approach (Juang et al., 2007;Luyssaert et al., 2014). The energy balance at the land surface is given by: where Rs i is incoming shortwave radiation (w m −2 ); Rl i is incoming longwave radiation (w m −2 ); Ɛ is emissivity; σ is Stefan-Boltzmann constant (W m −2 k −4 ); and Ts is LST. LE, H, and G are latent heat flux, sensible heat flux, and soil heat flux, respectively, as Eq. (3). Rearranging the formula and applying first-order derivative of the left-hand side and right-hand side of Eq. (9) gives the temperature change due to each of the energy balance terms: where changes in albedo (Δα), incoming shortwave radiation (ΔRs i ), incoming longwave radiation (ΔRl i ), latent heat flux (ΔλE), sensible heat flux (ΔH), ground heat flux (ΔG), and emissivity (ΔƐ) were computed using a 300 × 300 m grid for the corresponding land cover change pixels. For demonstration purposes, the energy flux changes and temperature decomposition were analyzed for Landsat 8 scene acquired on 2014-11-27.

Impacts of bushland loss on surface albedo, evapotranspiration, and surface temperature changes
The mean difference in albedo between adjacent cropland and bushland are displayed in Fig. 3a and b. The result shows a potential increase in albedo across all wavelength ranges (i.e., visible (VIS), near-infrared (NIR), and shortwave (SW)) if bushland is converted to cropland. The VIS albedo increases (P b 0.01) on average by 0.005 (SD 0.004) and 0.007 (SD 0.004) in Mwatate and Ndome, respectively. N95% of the investigated pixels indicated an average increase in VIS albedo. The average changes in NIR were non-significant (P = 0.46) and showed large variation in Mwatate (0.001, SD 0.01); in Ndome, changes in NIR were comparable to VIS despite the larger variation (0.007, SD 0.008).
The average SW broadband albedo change (P b 0.01) was 0.002 (SD 0.006) in Mwatate and 0.006 (SD 0.005) in Ndome, which was three times larger in magnitude. Consequently, the average surface radiative forcing (SRF) was stronger in Ndome (−6.02 W m −2 , SD 4.44) than in Mwatate (−2.38 W m −2 , SD 5.44). The average SRF showed great variation and can reach locally as high as −20 W m −2 ( Supplementary  Fig. S5).
Conversion of bushland to cropland leads to a potential increase in LST in both areas (P b 0.01; Fig. 3c). The average change was 1.75 K (SD 1.37) in Mwatate and 1.27 K (SD 1.12) in Ndome. From the total number of pixels affected by such conversion, 95% and 91% showed an increase in LST in Mwatate and Ndome, respectively. Furthermore, as the fraction of cropland and bushland varies from one pixel to another, the impact on LST difference was assessed (Fig. S6). In both areas, the increase or decrease in the cropland/bushland fraction had no impact (r 2 b 0.01, P N 0.05) on the LST difference ( Supplementary Fig. S6).
Bushland-to-cropland conversion causes a potential decrease in ET (P b 0.01; Fig. 3d). The decrease in ET is consistent with higher temperature and is built into the ET model. ET decreased on average by −0.45 mm day −1 (SD 0.29) in Mwatate and −0.41 mm day −1 (SD 0.31) in Ndome. The average change in ET was comparable in magnitude in the two areas. Locally, the maximum ET differences were as high as 1 mm day −1 in Mwatate and 0.7 mm day −1 in Ndome. The pixels that showed a decrease in ET accounted for 96% (in Mwatate) and 93% (in Ndome) of the total. Multiple-date testing also showed consistent decrease in ET during 2014 to 2018 ( Supplementary Fig. S7c, d).
In addition, the spatial, latitudinal, and longitudinal variation of ΔSW, ΔTs, and ΔET at 300 × 300 m resolution was displayed in Supplementary Fig. S9. The result showed that the local change was bigger in magnitude in all variables.
To understand the underlying physical mechanism causing the LST increase during bushland to cropland transitions, the changes in the surface energy balance components were assessed (Fig. 4). The result shows that the conversion of bushland to cropland tends to increase the upwelling shortwave radiation (SW u ) due to the higher albedo of cropland, increases the sensible and ground heat flux, and strongly reduces the latent heat flux (LE) in both areas (Fig. 4a, b). Given the arid water-limited characteristics of this area, these changes lead to elevated LST and ultimately increase the upwelling or emitted longwave radiation (LWu).

Relation of land surface temperature change with evapotranspiration and albedo changes
The relation of ΔLST with ΔET and albedo change (Δα) following cropland-to-bushland conversion was displayed in Fig. 5 using partial correlation and multiple regression. The ΔLST had a strong and negative correlation with ΔET (r~−0.8 maximum) in both areas (Fig. 5a) and the impact of ΔET on ΔLST was up to five times stronger than the Δα effect on ΔLST (Fig. 5b). The albedo change (Δα) had a weak and negative relation with ΔLST (r~−0.3 maximum). In general, ΔET and Δα together explained 59% and 70% of the variation in ΔLST in Mwatate and Ndome, respectively.
Since ΔET and Δα alone cannot fully explain the ΔLST, other factors that can affect ΔLST were also identified from the energy balance decomposition (see Section 2.3.4; Eq. (10)) and the contribution from each factor to ΔLST was quantified (Fig. 6). Both areas displayed similar patterns in all terms and the largest contribution to the total ΔLST came from the latent heat flux change (ΔLE) followed by sensible heat flux change (ΔH), ground heat flux change (ΔG), and upwelling longwave radiation flux change (ΔLW) (Fig. 6a, b). The contributions from albedo change (Δα), incoming shortwave radiation change (ΔSW), and emissivity change (ΔEM) were minor.
The net temperature change obtained from the sum of the individual components (modeled ΔLST, Eq. (10)) was compared with the observed ΔLST (from Landsat 8) to assess the robustness of the temperature decomposition (Fig. 7). The results show that the modeled ΔLST has good agreement with the observed ΔLST (r 2 = 0.95 for Mwatate and r 2 = 0.89 for Ndome).  Furthermore, the strong influence of ΔET on ΔLST was consistent across time, as demonstrated by the variation of LST and ET differences between bushland and cropland from multiple dates between 2014 and 2018 (Fig. 8). During this period, the ΔET explained 20% to 74% of the observed variation in ΔLST. The average ΔLST increases while the average ΔET decreases for all analyzed dates.

Discussion
Studies aiming to understand the impacts of land cover changes on regional and global climate have largely focused on forest ecosystems. On the other hand, the impacts of bushland clearing have thus far been poorly understood. Our results demonstrate that cropland expansion at the expense of bushland causes considerable local warming and at a magnitude comparable to the widely studied impacts of forest loss in the tropics. In the Horn of Africa, forest loss was shown to increase the surface temperature by up to 1.2 K, while locally the warming can be as high as 1.8 K (Abera et al., 2018). In comparison, our results reveal that bushland-to-cropland conversion can lead to an average warming of up to 1.75 K.
The clearing of bushland can have an important environmental impact at the regional level, given that a large part of the Horn of Africa is covered by Acacia-Commiphora bushland and thicket ecoregion (Bouvet et al., 2018;Dinerstein et al., 2017;White, 1983) and that the bushland is rapidly declining due to agricultural expansion, overgrazing, and fire Tadesse et al., 2014). The environmental impact of forest loss is relatively well known, and afforestation, reforestation, and conservation practices are consequently underway in the region as a part of REDD+ (Reducing emissions from deforestation and forest degradation and the role of conservation, sustainable management of forests and enhancement of forest carbon stocks in developing countries) (Mwangi et al., 2018). Furthermore, according to the highresolution (25 m) biomass map of African savannas and woodlands, there is quite a large amount of above-ground biomass (AGB) in this non-forest area (Bouvet et al., 2018). Therefore, AGB outside forest (including woodland and bushland) has been given more emphasis already from a carbon stock perspective. However, the conservation of Acacia-Commiphora bushland also requires attention considering their important impact on the LST. Therefore, with increasing population pressure on land, current and future clearing of bushlands for agriculture could exacerbate the warming and might hinder climate change mitigation efforts in the region. Below, we discuss the main biophysical mechanisms leading to temperature changes following bushlands clearing.

Impacts on albedo
Bushland-to-cropland conversion leads to an average increase in the broadband albedo in all spectral ranges (VIS, NIR, and SW). Multi-date testing of the broadband SW albedo changes during the period 2014 to 2018 also displayed a similar effect of such conversion (Supplementary Fig. S2c). This indicates that cropland has a higher albedo than  bushland on average. The lower albedo in bushlands, compared with cropland, can be explained by its darker appearance associated with its lower bare soil fraction, high canopy cover (40%), and greener canopy across different seasons (White, 1983; Supplementary Fig. S8). Unlike croplands, which mostly maintain bare soil cover during the dry seasons and sustain vegetation with shallower rooting systems, bushlands consist of plants (e.g., Acacia and Commiphora species) with deeper penetrating roots tolerant to arid conditions and appear greener and darker for longer periods than croplands.
When the spectral changes in albedo were compared, the changes in the VIS range were more consistent (i.e. N95% of pixels displayed an increase) and similar in magnitude than in the NIR and SW range. The  increase in VIS albedo after bushland-to-cropland conversion suggests reduced absorption of incoming radiation in the blue and red spectral regions associated with the decline in LAI and NDVI (Lillesand et al., 2008;Gibson, 2000;Supplementary Fig. S8). However, vegetation greening and LAI were usually not associated with a positive ΔNIR given the lower NDVI and LAI of croplands compared with bushlands ( Supplementary Fig. S8). In-situ studies in the region have shown that albedo intra-daily variability is dominated by changes in soil moisture. Thus, ΔNIR can be related to the background environmental conditions (i.e., soil color and moisture, amount of litter, and standing dead phytomass) as they are important in determining the albedo of an open canopy ecosystem (Samain et al., 2008;Campbell and Norman, 1998).
Even though the biophysical impacts of bushland-to-cropland conversion are relatively poorly studied, the effects of forest-to-cropland and other conversion classes were investigated in the region and beyond at coarser resolution (e.g., Abera et al., 2019;Duveiller et al., 2018;Bright et al., 2017). Similar to bushland-to-cropland conversions, these studies reported an average increase in SW albedo due to forestto-cropland conversion despite the structural and functional differences between forest and bushland. When the magnitude of albedo increase was compared, the bushland-to-cropland conversion had a smaller impact on albedo (maximum 0.6% increase) than forest-to-cropland conversion (1% increase) (Abera et al., 2019). In general, our result agrees with observational studies in the tropics that have reported surface brightening and consequent increase in albedo due to agricultural expansion at the expense of woody vegetation or taller to shorter vegetation transitions (Duveiller et al., 2018;Abera et al., 2019).

Impacts on evapotranspiration
The bushlands have higher ET rate than croplands; hence the conversion of bushland to cropland decreases ET. The higher ET rate in the bushlands can be related to the larger canopy cover (up to 40%), higher LAI (Fig. S8a, b), and relatively deep rooting system. Furthermore, the higher NDVI and LAI in the bushlands across 2014 to 2018 suggests higher photosynthetic and ET rates compared to croplands (Supplementary Figs. S7 and S8). Although the ET rate represents 24hour ET here, previous study showed that it is dominated by daytime ET as night time ET constitute smaller percentage (i.e. 1.7% of the 24hour ET; see the Malek (1992) for details).
The ΔET displayed a consistent decrease across multiple dates that were spread across the dry (e.g., January and February) and wet (e.g., May and November) seasons between 2014 and 2018. Despite the consistent decline in ET, however, a daily value does not represent monthly or seasonal variation in ΔET. Consequently, further studies are required at monthly and seasonal scales to determine if the decrease in ET due to bushland-to-cropland conversion persists across time scales. Future availability of eddy flux tower in the region can help to verify the magnitude of ΔET decline as well as inter-and intra -annual variabilities.

Impacts on surface energy balance and land surface temperature
The conversion of bushland to cropland perturbed the surface energy balance and resulted in a consistent increase in LST. The energy balance components responded differently to this perturbation in terms of magnitude of change. Although the increase in reflected shortwave radiation (due to an albedo increase) causes a local cooling effect, the reduction of latent heat flux was comparatively stronger and outweighed the albedo effect by approximately five times on average, thus leading to net warming (Fig. 6a, b). This consequently led to the increase in emitted longwave radiation and sensible heat flux. Similar to our result, several global and regional observational and modeling studies reported the dominant influence of ET on LST in the tropics during conversion of woody vegetation, such as forest to cropland or grassland (e.g., Abera et al., 2019;Duveiller et al., 2018;Bright et al., 2017;Alkama and Cescatti, 2016;Davin and de Noblet-Ducoudre, 2010).
For instance, the albedo cooling effect due to conversion of forest to cropland in the Horn of Africa was outweighed by non-radiative processes (ET and surface roughness) by up to 10 times (Abera et al., 2019). Similarly, a global study showed that the warming due to vegetation changes in the tropics was mainly associated with agricultural expansion, and the increase in temperature due to the reduction of ET was not counterbalanced by the increase in albedo (surface brightening) (Duveiller et al., 2018).
The mechanism behind the reduction in ET and the consequent warming can be related to the ET and surface roughness changes associated with the reduction of woody cover during bushland-to-cropland transitions (Rotenberg and Yakir, 2010;Pitman, 2003;Hoffmann and Jackson, 2000;Xue et al., 1996). The decline in woody cover will reduce ET rate as indicated across 2014 to 2018 ( Supplementary Fig. S8c, d). Compared to croplands, trees (e.g., Acacia and Commiphora) have better efficiency in transferring water from deep soil to the atmosphere through their larger leaf area and deeper rooting systems (Davin and de Noblet-Ducoudre, 2010).

Spatio-temporal biophysical variability
The biophysical changes (ΔSW broadband albedo, ΔET, and ΔLST) as well as the relationship between ΔET and ΔLST showed large spatiotemporal variations during bushland-to-cropland conversion. The spatial variability is likely related to the background environmental conditions (e.g. soil moisture variabilities) and bare soil fraction (Abera et al., 2020). The stronger soil-moisture variability in semi-arid environments (Jung et al., 2010) can affect the magnitude of ΔSW, ΔET, and ΔLST across the study area (Abera et al., 2020).
The inter-annual variabilities (Figs. S2c and S7) are likely related to global and regional climate change factors. For instance, associated with the recurrent rainfall extremes in the region, the SW albedo can increase during drought events whereas the opposite pattern happens during extreme wet periods (Abera et al., 2020). This also affects the surface-atmosphere energy exchange and hence ET and LST (Teuling et al., 2010). Increase in CO2 concentration may affect ET rate through their impact on water use efficiency (Hatfield and Dold, 2019;Eamus, 1991). These factors in turn can influence the ΔET and ΔLST relationship through impacting transpiration and energy storage in soil and plants (Sun et al., 2016).
Impacts of local anthropogenic factors can be minimum as stable adjacent land cover pixels were investigated and agriculture mainly depends on seasonal rainfall in the study area. For full understanding of the causes of the inter-annual variabilities of ΔET-ΔLST relationship, future study can benefit from exploring various atmospheric and surface biophysical factors as well as their interactions comprehensively.

Limitations and uncertainties
The main limitations and uncertainties of our assessment are associated with lack of data availability and validation, and model requirements. The biophysical changes (ET and albedo) and consequent impact on LST were assessed on a daily time scale corresponding to the cloud free Landsat 8 satellite images. The changes on a monthly, seasonal, and annual time scale were not assessed due to the limitations of the Landsat 8 temporal resolution and persistent cloud cover in the region.
Comparison of ΔLST from Landsat 8 and modeled ΔTs from energy balance decomposition (EBD) needs to be independent. In the presented form, LST from Landsat 8 was also applied in the EBD (Eq. (10)) and hence our EBD results can carry uncertainties. For this, future availability of upwelling longwave radiation measurement from flux tower can help to evaluate the performance of EBD approach in the region. However, Eq. (10) contains energy balance components which were derived from Landsat 8 bands other than thermal band, weather station, and digital elevation data in METRIC model. In that sense, the result presented in our study can be used as an indicator to measure the performance of the EBD approach in areas of scarce ground flux data. Although validation of the ET estimation from METRIC modeling was not possible due to lack of ground flux data in this region, intermodel comparison was made between METRIC model output and Penman-Monteith model output which was computed using the Voi station data during 2014-2018 (Supplementary Table S2). The point comparison at Voi station showed a mean difference (MD) of −0.05 mm day −1 and mean absolute difference of 0.47 mm day −1 . This difference is within the range of previous estimates in the region using empirical models (Maeda et al., 2011).
Therefore, further studies are needed at monthly, seasonal, and annual time-scales to get complete understanding of the biophysical changes and their associated impact on the climate following bushland to cropland conversion. In this regard, availability of flux tower data in the future can give better quantification of the actual impacts and validation of satellite-based studies.

Conclusions
We estimated the potential impacts of bushland-to-cropland conversion on surface biophysical properties (albedo and ET) and the consequent effect on local climate through LST changes in two study sites in southern Kenya. The conversion leads to a consistent increase in spectral albedo but the change in the shortwave range was dominated by the change in the visible range. The increase in albedo reduces the available shortwave energy at the surface by up to 6.02 W m −2 and generates an average cooling of −0.8 K. Nevertheless, such cooling was shown to be outweighed by ET decline, resulting in an average warming of up to 1.75 K. This warming due to bushland-to-cropland conversion has considerable impact on the environment due to its comparable magnitude with the effects of forest loss from previous study (Abera et al., 2018) in the region.
This finding reveals the potential of Acacia-Commiphora bushlands in regulating LST in the region and suggests the necessity of their conservation. Furthermore, it highlights that clearing of this vegetation for cropland induces warming. Hence, current and future clearing of bushlands for agriculture need to carefully consider their climatic impact.

Declaration of competing interest
The authors certify that they have NO affiliations with or involvement in any organization or entity with any financial interest (such as honoraria; educational grants; participation in speakers' bureaus; membership, employment, consultancies, stock ownership, or other equity interest; and expert testimony or patent-licensing arrangements), or non-financial interest (such as personal or professional relationships, affiliations, knowledge or beliefs) in the subject matter or materials discussed in this manuscript.