Expansion of oil palm and other cash crops causes an increase of land surface temperature 1 in Indonesia 2 3

23 24 Indonesia is currently one of the regions with the highest transformation rate of the land surface 25 worldwide due to the expansion of oil palm plantations and other cash crops replacing forests 26 on large scales. Land cover changes, which modify land surface properties, have a direct effect 27 Biogeosciences Discuss., doi:10.5194/bg-2017-203, 2017 Manuscript under review for journal Biogeosciences Discussion started: 29 May 2017 c © Author(s) 2017. CC-BY 3.0 License.

S1. Surface temperature retrieval from Landsat thermal images Table S1.1.Steps in the retrieval of the surface temperature from Landsat TIR band Table S1.2.LMIN and LMAX values for Landsat 7 ETM+ Table S1.3.Mean solar exo-atmospheric irradiance (ESUNλ) for Landsat 7 ETM+ S2.Atmospheric correction of the thermal band Table S2.1.Input and output parameters for/from NASA´s online atmospheric correction parameter calculator S3.ET from satellite images with SEBAL Fig. S3.1 Analysis of the steps involved in deriving the input for deriving ET from Landsat images with SEBAL Fig. S3.2Comparison of ET derived from upper anchor and lower anchor pixels.Table S3.1.u*, rah, LE and H measured at a young and mature oil palm plantation S4.Mean LST, NDVI, Albedo and NDVI extracted for 7 land cover types Fig. S4.1 Mean LST, NDVI, Albedo and NDVI extracted from Landsat LST images for 7 land cover types S5.Difference in LST, NDVI, albedo and ET between Forest (FO) and 6 other land cover types Fig. S5.1 Differences in LST (∆LST), NDVI (∆NDVI), Albedo (∆Albedo) and Evapotranspiration (∆ET) between other land covers (RU, MOP, PF, YOP, UB and CLC) and forest (FO) in the Jambi province S6.Statistical analysis Table S6.1 ANOVA statistics Table S6.2Post-hoc Tukey HSD test statistics Table S6.S9.1 Land use change (1990) -2000 -2010 Table S9.2Contribution of land cover change to total LST increase S1.Surface temperature retrieval from Landsat thermal images Surface temperature from Landsat ETM+ is derived in a series in steps using the red (R), near-infrared (NIR) and thermal infrared (TIR) band, and follows the method described by Bastiaanssen et al. (1998aBastiaanssen et al. ( , 1998b) ) and summarized in table S1.  sensor constants for converting band 6 to surface temperature mW/cm 2 /sr/μm K #: when using Landsat surface reflectance product, this step is not necessary.§The transformed SAVI is considered here to be equal to the LAI.References: 4: Huete (1988); 5: Bulcock and Jewitt (2010);6, 7: Bastiaanssen et al. (1998a6, 7: Bastiaanssen et al. ( , 1998b)); 8: Wukelic et al. (1989); Coll et al. (2010).

S2. Atmospheric correction of the thermal band
For the atmospheric correction of the thermal radiance, the correction parameters τ, Rp and Rsky are derived from NASA's online atmospheric correction calculator (Barsi et al., 2003(Barsi et al., , 2005) ) which requires the following input: latitude and longitude, elevation, air temperature, surface pressure and relative air humidity.This Web-based ACT has been developed for TM and ETM+ thermal data (Barsi et al., 2005;Coll et al., 2010).It uses atmospheric profiles from the National Centers for Environmental Prediction (NCEP) interpolated to a particular location, date and time, and the MODTRAN-4 code to calculate the atmospheric-correction parameters for the bandpass of either the TM or ETM+ thermal band for a given date and site (Coll et al., 2010).The output of the online atmospheric parameter calculator is: band average atmospheric transmission (τ), effective bandpass upwelling radiance (Rp) and the effective bandpass downwelling radiance (Rsky).The calculated atmospheric parameters can be applied to a given scene to retrieve the surface temperature for the area of interest.The thermal band (L6) is corrected for atmospheric effects after Wukelic et al. (1989) and Coll et al. (2010) as: The input data required for the atmospheric calculator and the output parameters required for atmospheric correction and applied to the selected Landsat satellite image are summarized in table S2.1.4.25 Band average atmospheric transmission = τ, effective bandpass upwelling radiance = Rp (W/m 2 /sr/µm), effective bandpass downwelling radiance = Rsky (W/m 2 /sr/µm).
εa: the atmospheric emissivity (-), estimated with: εa = 1 -0.26 • exp {-7.77 × 10 -4 ) • (273.15-Ta) 2 } (Eq.S2) σ: Stephan-Boltzmann constant (5.67 × 10 -8 W/m 2 /K 4 ); LST: the surface temperature (K) derived from Landsat; Ta: is the near surface air temperature (K).We calculated u* and rah for a young and mature oil palm plantation (table S3.1).These were calculated from meteorological measurement on these locations.Because the satellite image was acquired outside the time period in which meteorological measurements were made, we selected dates and times that had similar conditions on the day the image was acquired.We used incoming shortwave radiation as a main criteria (> 690 Wm -2 and < 720 Wm -2 , between 10.00 and 11:00 am local time).u* and rah derived from the satellite image show a certain level of agreement with the u* and rah calculated from meteorological data.
4. Latent heat flux (LE, W/m 2 ) is estimated as residual from Net radiation, Ground heat and sensible heat flux as: LE = Rn -G -H (Eq.S11) 5. Instantaneous evapotranspiration (ET, mm/hr) for each pixel is estimated from LE as: = 3600 (Eq.S12) 3600 is the time conversion from seconds to hours λ = latent heat of vaporization (2.43×10 6 J/kg) We tested different combinations of 5 hot and 5 cold pixels, that could serve as anchor pixels, and then compared the effects of anchor pixel selection on the ET output.Our comparison showed that the anchor pixels we selected showed an overall effect on the magnitude of ET of less than 10% and had no effect on the ranking of the ET by land use type (Fig. S3.2).
We calculated LE and H for a young and mature oil palm plantation (table S3.1).These were calculated from flux measurements on these locations.Because the satellite image was acquired outside the time period in which meteorological measurements were made, we selected dates and times that had similar conditions on the day the image was acquired.LE and H derived from the satellite image show some agreement with the LE and H calculated from meteorological data.
Table S3.H.These values were calculated for the lower (> 690 W/m 2 ) and upper (< 720 W/m 2 )limit of incoming solar radiation we selected to match conditions on the day the satellite image was acquired.Both in situ and MODIS observations are consistent, i.e. the morning temperatures (10:30 am local time vs Terra LST) were warmer than the evening temperatures (10:30 pm local time vs Aqua LST).Differences between the two sources are caused by the comparison of point measurements with pixel values, differences in spatial resolution, differences in soil contribution to the LST estimate, distance in LST measurements and particularly differences in emissivity used for temperature correction.The thermal infrared sensor measuring the surface canopy temperature of the oil palm plantation had fixed default values.MODIS emissivity is derived from 3 thermal bands and adjusted accordingly for every measurement.Remarks/explanation: *) units in fraction of the land cover, which is calculated as the observed change (%, table S9.1) divided by 100 **) warming effect from Fig. 4a (main text) 1) From the observed changes, but seems to be very underestimated, we assume the remaining rubber is mixed with forests as Jungle Rubber 2) MOP = we assume half of the plantations to be mature 3) Not included because increase in area is not known 4) YOP = we assume half of the plantations to be young 5) We assume half of "other" to be urban areas (this assumption is already an overestimation) 6) We assume 50% of the bush to be degraded areas and barren and 50% of the class other to be clear cut land 7) 50% of class "bush" is not included, because this class was not part of our study

3
The relation LST-Albedo-NDVI-ET separated by land cover type S7.Comparison of MODIS LST to in situ measured canopy LST Fig. S7.1 MODIS LST compared to in situ measured canopy surface temperature.S8.Comparison of MODIS Air temperature with locally measured air temperature Fig. S8.1 MODIS Air temperature compared with in situ measured air temperatures S9.Land use change analysis for the Jambi province for 2000 -2010 Table Fig. S3.1 Analysis of the steps involved in deriving the input for deriving ET from Landsat images with SEBAL.
Fig. S4.1 Mean LST, NDVI, Albedo and NDVI extracted from Landsat LST images for 7 land cover types.The values were extracted from small plots that could only be used for Landsat images.
Fig. S7.1 MODIS LST compared to in situ measured canopy surface temperature.Canopy surface temperature was measured above a homogeneous mature oil palm plantation (12 years old).MODIS LST during 1½ year (mid 2014 -end 2015) was extracted from the pixel covering the location where the in situ canopy surface temperature was measured.LST from MODIS on the Terra and Aqua platform were used: Aqua LST in our comparison were measured in the evening hours (around 22:30, local time), Terra LST were measured during morning hours (10:30 am, local time).

Fig
Fig. S8.1 MODIS Air temperature compared with in situ measured air temperatures Young oil palm plantation (2 years), Urban area and Open field (surrounded by forest) are land use types where the meteorological towers in Jambi were located.MODIS air temperatures were extracted from the MODIS Air temperature profile product (MOD07) for the three locations.

Table S1 .
1. Steps in the retrieval of the surface temperature from Landsat 7 TIR band

Table S2
.1.Input and output parameters for/from NASA´s online atmospheric correction parameter calculator applied to two satellite images of the study area.Results in this study were based on the satellite image acquired on 2013-06-19.