Expansion of oil palm and other cash crops causes an increase of the land surface temperature in the Jambi province in Indonesia

Indonesia is currently one of the regions with the highest transformation rate of land surface worldwide related to the expansion of oil palm plantations and other cash crops replacing forests on large scales. Land cover changes, which modify land surface properties, have a direct effect on the land surface temperature (LST), a key driver for many ecological functions. Despite the large historic land transformation in Indonesia toward oil palm and other cash crops and governmental plans for future expansion, this is the first study so far to quantify the impacts of land transformation on the LST in Indonesia. We analyze LST from the thermal band of a Landsat image and produce a highresolution surface temperature map (30 m) for the lowlands of the Jambi province in Sumatra (Indonesia), a region which suffered large land transformation towards oil palm and other cash crops over the past decades. The comparison of LST, albedo, normalized differenced vegetation index (NDVI) and evapotranspiration (ET) between seven different land cover types (forest, urban areas, clear-cut land, young and mature oil palm plantations, acacia and rubber plantations) shows that forests have lower surface temperatures than the other land cover types, indicating a local warming effect after forest conversion. LST differences were up to 10.1± 2.6 C (mean±SD) between forest and clear-cut land. The differences in surface temperatures are explained by an evaporative cooling effect, which offsets the albedo warming effect. Our analysis of the LST trend of the past 16 years based on MODIS data shows that the average daytime surface temperature in the Jambi province increased by 1.05 C, which followed the trend of observed land cover changes and exceeded the effects of climate warming. This study provides evidence that the expansion of oil palm plantations and other cash crops leads to changes in biophysical variables, warming the land surface and thus enhancing the increase of the air temperature because of climate change.


Introduction
Indonesia is one of the regions where the expansion of cash crop monocultures such as acacia (timber plantations), rubber, oil palm plantations and smallholder agriculture has drastically reduced the area of primary forest in the last 2.5 decades (Bridhikitti and Overcamp, 2012;Drescher et al., 2016;Marlier et al., 2015;Miettinen et al., 2012;Verstraeten et al., 2005).This large-scale conversion of rainforest for agricultural use has been observed on the island of Sumatra, which has experienced the highest primary rainforest cover loss in all of Indonesia (Drescher et al., 2016;Margono et al., 2012;Miettinen et al., 2011).Forest cover in the Sumatran provinces of Riau, North Sumatra and Jambi declined from 93 to 38 % of provincial area between 1977 and2009 (Mietti-Published by Copernicus Publications on behalf of the European Geosciences Union. nen et al., 2012).These large-scale transformations, observed as land cover change, and land use intensification have led to substantial losses in animal and plant diversity, ecosystem functions and changed microclimatic conditions (Clough et al., 2016;Dislich et al., 2016;Drescher et al., 2016).Additionally, these changes directly alter vegetation cover and structure and land surface properties such as albedo, emissivity and surface roughness, which affect gas and energy exchange processes between the land surface and the atmosphere (Bright et al., 2015).
Replacing natural vegetation with another land cover modifies the surface albedo, which affects the amount of solar radiation that is absorbed or reflected and consequently alters net radiation and local surface energy balance.A lower or higher albedo results in a smaller or greater reflection of shortwave radiation.As a result, the higher or lower amounts of net radiation absorption may increase or decrease the surface temperature and change evapotranspiration (ET) (Mahmood et al., 2014).
Changes in land cover also alter surface emissivity, i.e. the ratio of radiation emitted from a surface to the radiation emitted from an ideal black body at the same temperature following the Stefan-Boltzmann law.Emissivity of vegetated surfaces varies with plant species, density, growth stage, water content and surface roughness (Snyder et al., 1998;Weng et al., 2004).A change of emissivity affects the net radiation because it determines the emission of longwave radiation that contributes to radiative cooling (Mahmood et al., 2014).
Water availability, surface type, soil humidity, local atmospheric and surface conditions affect the energy partitioning into latent (LE), sensible (H) and ground heat (G) fluxes (Mildrexler et al., 2011).Surface roughness affects the transferred sensible and latent heat by regulating vertical mixing of air in the surface layer (van Leeuwen et al., 2011), thereby regulating land surface temperature (LST).Through its association with microclimate, net radiation and energy exchange (Coll et al., 2009;Sobrino et al., 2006;Voogt and Oke, 1998;Weng, 2009;Zhou and Wang, 2011), LST is a major land surface parameter, and as a climatic factor it is regarded to be a main driver of diversity gradients related to the positive relationships between temperature and species richness (Wang et al., 2016).
The replacement of natural vegetation also changes ET (Boisier et al., 2014) and LST because the surface biophysical variables (i.e.surface albedo, LST, emissivity and indirectly leaf area index (LAI) and normalized difference vegetation index (NDVI)) are interconnected through the surface radiation balance.When ET decreases, for example, surface temperatures and sensible heat fluxes increase; in contrast, when ET increases, the increased LE fluxes lower surface temperatures and decrease H fluxes (Mahmood et al., 2014) under equal net radiation conditions because with a change in vegetation, soil and ecosystem heat flux and net radiation also change due to an alteration of the biophysical variables.Vegetation structure, represented by NDVI, LAI and vegeta-tion height, is in this respect an important determinant of the resistances or conductivities to heat, moisture and momentum transfer between the canopy and the atmosphere (Bright et al., 2015), facilitating the amounts/ratios of sensible heat to water vapor dissipation away from the surface (Hoffmann and Jackson, 2000).
To understand the effects of land cover changes on LST, the associated biophysical variables must be evaluated.This can be done through the surface radiation budget and energy partitioning which unite these biophysical variables directly or indirectly: albedo as direct determinant of the net solar radiation, NDVI as a vegetation parameter determining the emissivity, which in turn determines the amount of reflected and emitted longwave radiation; LST directly affecting the amount of emitted longwave radiation from the surface; and ET, which affects the amount of energy that is used for surface cooling via the evaporation of water.
The effect of land cover change on LST is dependent on the scale, location, direction and type of the change (Longobardi et al., 2016).Several studies showed an LST increase after forest conversion to built-up areas, agricultural land (Zhou and Wang, 2011), crop land and pasture lands (Peng et al., 2014) in China.Similar observations were reported for South American ecosystems: low vegetation such as grasslands in Argentina were warmer than tall tree vegetation (Nosetto et al., 2005).In Brazil, the surface temperature increased after the conversion of natural Cerrado vegetation (a savanna ecosystem) into crop/pasture (Loarie et al., 2011a).Similar effects were also observed for other South American biomes (Salazar et al., 2016).In a global analysis, Li et al. (2015) showed that the cooling of forests is moderate at midlatitudes but northern boreal forests are warmer, an indication that the effect of land cover change on LST varies with the location of the land cover change (Longobardi et al., 2016).Similar studies on the Indonesian islands are lacking but surface temperature increases are expected as an effect of oil palm and cash crop land expansion in the recent decades.
Measuring LST changes is critical for understanding the effects of land cover changes, but challenging.LST can be monitored with LST products retrieved from thermal infrared (TIR) remote sensing data: e.g. the use of the thermal bands of the Moderate Resolution Imaging Spectrometer (MODIS) on board the Terra and Aqua satellite (Sobrino et al., 2008), the thermal band of the Thematic Mapper (TM) on board the LANDSAT-5 platform (Sobrino et al., 2004(Sobrino et al., , 2008) ) or Enhanced Thematic Mapper (ETM+) on board the LANDSAT-7 platform.The advantage of MODIS data is the availability of readily processed products at high temporal resolution (daily) at medium (250-500 m) to coarse (1000-5000 m) spatial resolution scale; MODIS LST product (MOD11A1/MYD11A1), for example, is provided at a daily temporal resolution with a spatial resolution of 1 km.Landsat data are provided at a higher spatial resolution (30 m), but the temporal resolution is limited to 16 days and the retrieval of LST requires the correction of the satellite-observed radiances for atmospheric absorption and emission (Coll et al., 2009).Besides LST, the connected biophysical variables of the energy and radiation budget can be derived from the visible and near-infrared (VIS-NIR) bands of MODIS or Landsat, making integrated monitoring of the biophysical variables related to changing land surface possible.In Indonesia, a large proportion of the land use changes is driven by smallholders (Dislich et al., 2016), and thus a combination of Landsat (for a fine spatial resolution) and MODIS (for temporal developments) seems desirable.
The modification of the physical land surface properties influences climate and local microclimatic conditions via biogeochemical and biophysical processes.Therefore, given Indonesia's history of large-scale agricultural land conversion and governmental plans to substantially expand the oil palm production (Wicke et al., 2011), it is important to study the effects of the expansion of cash crop areas on the biophysical environment, especially on LST as a key land surface parameter.These effects have been poorly studied in this region and, according to our knowledge, this is the first study to quantify the effects of land use change on LST in Indonesia.We focus on the Jambi province (in Sumatra, Indonesia) as it experienced large land transformation towards oil palm and other cash crops such as rubber plantations in the past, and it may serve as an example of future changes in other regions.
Our main objective is to quantify the differences in LST across different land cover types and to assess the impact of cash crop expansion on the surface temperature in the Jambi province in the past decades (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015).With this study we aim (1) to evaluate the use of Landsat and MODIS satellite data as sources of a reliable surface temperature estimation in a tropical region with limited satellite data coverage by comparing the surface temperatures retrieved from both satellite sources to each other and against ground observations, (2) to quantify the LST variability across different land cover types, (3) to assess the long-term effects of land transformation on the surface temperature against the background of climatic changes and (4) to identify the mechanisms that explain the surface temperature changes caused by alterations of biophysical variables.In this study we compare the surface temperatures of different land cover types that replace forests (i.e. oil palm, rubber and acacia plantations, clear-cut land and urban areas) by using high-resolution Landsat and medium-resolution MODIS satellite data and discuss the differences by taking into account other biophysical variables such as the albedo, NDVI and ET.

Study area
The study was carried out in the lowlands (approx.25 000 km 2 ) of the Jambi province (total area 50 160 km 2 ) on Sumatra, Indonesia, between latitudes 0 • 30 and 2 • 30 S and longitudes 101 • and 104 • 30 E (Fig. 1).This region has undergone large land transformation towards oil palm and rubber plantations over the past decades and thus may serve as an example of expected changes in other regions of Indonesia (Drescher et al., 2016).The area has a humid tropical climate with a mean annual temperature of 26.7 ± 0.2 • C (1991-2011, annual mean ± SD of the annual mean), with little intra-annual variation.Mean annual precipitation was 2235 ± 381 mm and a dry season with less than 120 mm monthly precipitation usually occurred between June and September (Drescher et al., 2016).Previously logged rainforests in the Jambi province have been converted to intensively managed agro-industrial production zones and into smallholder farms to grow cash crops of rubber (Hevea brasiliensis) and oil palm (Elaeis guineensis) or fast-growing tree species such as Acacia mangium for pulp production (Drescher et al., 2016).The area cultivated with oil palm grew faster than the area cultivated with rubber plantations between 1990 and 2011 (Clough et al., 2016).
For this study, we used two data sets of different plot sizes.For the first data set, we delineated 28 large plots (ranging from 4 to 84 km 2 ) of 7 different land cover types: forest (FO), rubber (RU), acacia plantation forest (PF), young oil palm plantation (YOP), mature oil palm plantation (MOP), urban area (UB) and clear-cut areas (CLC) (Fig. 1).The delineation was based on visual interpretation in combination with field observations, which were conducted between October and December 2013.The large size of the plots was necessary to make a comparison between MODIS and Landsat images (see Sect. 2.3).For the second data set, we selected 49 smaller plots within and outside these 28 large plots (between 50 m × 50 m and 1000 m × 1000 m) (Fig. 1), which allowed us to increase the number of plots to use when analyzing Landsat images.These small plots were used to extract the LST, NDVI, albedo (α) and ET from a high-resolution Landsat satellite image (see section satellite data) for the which different land cover types of interest.

Meteorological data
Air temperature and relative air humidity were measured at four reference meteorological stations located in open areas within the study area (Drescher et al., 2016), with thermohygrometers (type 1.1025.55.000,Thies Clima, Göttingen, Germany) placed at 2 m height.Measurements were recorded every 15 s and then averaged and stored in a DL16 Pro data logger (Thies Clima, Göttingen, Germany) as 10 min mean, from February 2013 to December 2015.We used the air temperature from the meteorological stations to compare to MODIS air temperatures (MOD07_L2).The relative air humidity was used as an input parameter for NASA's online atmospheric correction (ATCOR) parameter tool to derive parameters to correct Landsat thermal band for atmospheric effects (see Sect. 2.3).We also used air temperature and relative humidity (RH) from two eddy covariance flux towers located in the study area (Meijide et al., 2017): one in a young oil palm plantation (2 years old; 01 • 50.127 S, 103 • 17.737 E) and the other one in a mature oil palm plantation (12 years old; 01 • 41.584 S, 103 • 23.484 E).At these flux towers, air temperature and relative humidity were measured above the canopy with the same instruments as in the reference meteorological stations (see Meijide et al., 2017, for a description of the methodology).At the flux tower located in the mature oil palm plantation, we also measured the surface canopy temperature between August 2014 and December 2015, which was compared to MODIS LST estimates from the same period.The canopy temperature was measured with two infrared sensors (IR100) connected to a data logger (CR3000), both from Campbell Scientific Inc. (Logan, USA).For a regional coverage we used ERA-Interim daily air temperature grids (http://apps.ecmwf.int/datasets/data/interim-full-daily/levtype=sfc/; Dee et al., 2011) from 2000 to 2015 at 0.125 • resolution to study the annual air temperature trend in this period.(1)

Land surface temperature
LST was derived following the method proposed by Bastiaanssen (2000), Bastiaanssen et al. (1998a), Coll et al. (2010) and Wukelic et al. (1989) for computing the surface temperature from the TIR band (band 6) of Landsat (Supplement S1).The TIR band was first converted to thermal radiance (L6, W m −2 sr −1 µm −1 ) and then to atmospherically corrected thermal radiance (R c , W m −2 sr −1 µm −1 ), as described by Wukelic et al. (1989) and Coll et al. (2010), and with the atmospheric parameters obtained on NASA's online Atmospheric Correction Calculator (Barsi et al., 2003(Barsi et al., , 2005) ) (Supplement S2).The LST (K) was computed with the following equation similar to the Planck equation, as in Coll et al. (2010) and Wukelic et al. (1989): where εNB is the emissivity of the surface obtained from the NDVI (Supplement Table S1), k1 (= 666.09 mW cm −2 sr −1 µm −1 ) and k2 (= 1282.71K) are sensor constants for converting the thermal radiance obtained from band 6 of Landsat 7 to surface temperature.The surface temperature derived from Landsat thermal band was compared with the MODIS LST product that was acquired on the same day at 10:30 local time.The Landsat LST image was first resampled to MODIS resolution to enable a pixel-to-pixel comparison, followed by extracting the average LST of 7 land cover types with the data set containing the large delineated plots (Fig. 1).

Evapotranspiration
Based on the Surface Energy Balance Algorithm for Land (SEBAL) (Bastiaanssen, 2000;Bastiaanssen et al., 1998a, b) we estimated ET (mm h −1 ) from latent heat fluxes (LE, W m −2 ), which were computed as the residual from sensible (W m −2 ) and ground (W m −2 ) heat fluxes subtracted from net radiation (R n , W m −2 ): We calculated R n as the sum of incoming shortwave and longwave radiation minus the reflected shortwave and longwave radiation and the emitted longwave radiation (Eq.5).
The surface albedo, surface emissivity and surface temperature determine the amounts of incoming and reflected radiation: where S d ↓ is the incoming shortwave solar radiation (W m −2 ) at the surface, α is the surface albedo (Eq.2); ε 0 is the surface emissivity (−), ε a is the atmospheric emissivity (−), σ is the Stefan-Boltzmann constant (5.67 × 10 −8 W m −2 K −4 ), LST is the surface temperature (K, Eq. 3) and T a is the sky temperature (K).The surface emissivity (ε 0 ) is derived from the NDVI and is described in the Supplement (Table S1).The average atmospheric emissivity (ε a ) is estimated with the model of Idso and Jackson (1969): Ground heat fluxes (W m −2 ) were derived as a fraction of R n from an empirical relationship between LST, α and NDVI (Bastiaanssen, 2000) as In SEBAL sensible heat (W m −2 ) was calculated as where ρ is the air density (1.16 kg m −3 ); C p is the specific heat of air at constant pressure (1004 J kg −1 K −1 ); r ah is the aerodynamic resistance to heat transport (s m −1 ); and a and b are regression coefficients, which are determined by a hot extreme pixel (where LE = 0 and H is maximum) and a cold extreme pixel (where H = 0 and LE is maximum).The aerodynamic resistance to heat transport, r ah , is calculated through an iterative process with air temperature measured at 2 m as input.SEBAL is described in Bastiaanssen (2000) and Bastiaanssen et al. (1998a, b).The application of SEBAL in this research is briefly described in the Supplement S3.

Local short-term differences between different land cover types
From the created LST, NDVI, albedo and ET images we extracted the average values of the different land cover classes with the data set containing the small 49 delineated plots covering 7 different land cover types (Fig. 1).The average effect of land transformation, i.e. the change from forest to another non-forest land cover type, on the surface temperature was evaluated as (see Li et al., 2015) A negative LST indicates a cooling effect and positive LST indicates a warming effect of the non-forest vegetation compared with forest.The same procedure was applied in evaluating the effect of land transformation on the NDVI, albedo and ET.

Effects of land cover change on the provincial surface temperature in the past decades
To analyze the long-term effects on the provincial scale we used the MODIS daily LST time series (MOD11A1 and MYD11A1) from 2000 to 2015.MOD11A1 provides LST for 10:30 and 22:30 and we used the times series between 2000 and 2015.MYD11A1 provides LST for 01:30 and 13:30 and is available from 8 July 2002; we used complete years in our analysis and therefore used the MYD11A1 time series from 2003 to 2015.We calculated the mean annual LST at four different times of the day (10:30, 13:30, 22:30 and 01:30) between 2000 and 2015 for the lowland of Jambi from the MODIS daily LST time series (MOD11A1 and MYD11A1).First, we calculated for each pixel the average LST pixel value using only the best-quality pixels for every year.Then, from these pixels we made a composite image (n = 16, one for each year) for the province.Finally, from each composite image we calculated the mean annual lowland provincial temperature as the average of all the pixels that are enclosed by a zone delineating the lowland of the Jambi province.We performed the same analysis with the MODIS 16-day NDVI product (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) and the ERA daily temperature grid (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) to compare the annual trends of LST, NDVI and air temperature of the province.The average provincial LST and NDVI were compared with the mean LST and NDVI of a selected forest that remained undisturbed forest during the 2000-2015 period.

Statistical analysis
For a comparison of the Landsat-derived LST and the MODIS LST we analyzed the statistical relationships with the coefficient of determination (R 2 ), the root mean square error (RMSE), the mean absolute error (MAE) and the bias: where O i is MODIS LST, E i is the Landsat surface temperature and N is the number of pixels compared.Model type 2 linear regression was applied for fitting the relation between MODIS LST and Landsat LST.
We tested the relation between the biophysical variables LST (or L6 and R c , both as pre-or intermediate products before obtaining LST), albedo (α), NDVI and ET with a correlation analysis and a multiple linear regression was applied to analyze the effects of the biophysical variables on the LST.We used the model LST (or R c or L6) ∼ α+ NDVI + ET and used R 2 and standardized β coefficients to evaluate the strength of the biophysical variables in predicting the LST.

Landsat LST compared to MODIS LST
Landsat and MODIS images showed similar spatial LST patterns (Fig. 2).In both images the relatively hot areas (red) correspond to the known clear-cut areas, urban areas or other sparsely vegetated areas, the relatively cool areas (blue) correspond to vegetated areas such as forest, plantation forests and mature oil palm plantations.The coarse-resolution scale of MODIS (1000 m for LST) allows a large regional coverage of the study area but does not allow to retrieve detailed information on small patches (smaller than 1 km 2 ).In contrast, the Landsat 7 image allows a detailed study of patches that are small enough (as small as 30 × 30 m 2 ) but it is affected by the scan line error causing data loss at the edges of the image.In both MODIS and Landsat images clouds and cloud shadows were removed and therefore lead to data gaps in the images.
Similar differences were found for the NDVI between forest and the other land covers (Fig. 4b).The negative NDVI indicates that the non-forest land cover types had lower NDVI than the forest.NDVI between FO and RU, MOP, PF and YOP were small (between −0.01 ± 0.02 ( NDVI MOP − FO ) and −0.12 ± 0.06 ( NDVI YOP − FO )).The largest NDVIs were between forest and the non-vegetated land cover types, i.e.UB and CLC ( NDVI = −0.42± 0.11 and −0.41 ± 0.08, respectively).
A separate analysis (Table S6.3, Supplement S6) showed that ET was a strong predictor of LST for each land cover type in this study and that NDVI and albedo were minor predictors of LST.

Effects of land use change on the provincial surface temperature in the past decades
The average annual LST of Jambi was characterized by a fluctuating but increasing trend during daytime (Fig. 5a  and b) between 2000 and 2015.The average morning LST (10:30) increased by 0.07 • C per year (R 2 = 0.59; p < 0.001), while the midday afternoon LST (13:30 local time) increased by 0.13 • C per year (R 2 = 0.35; p = 0.02) between 2003 and 2015.While the daytime LST showed a clear increase, the night and evening LST (22:30 and 01:30, Fig. 5c and d) trends showed a small decrease of −0.02 • C (R 2 = 0.29; p = 0.02) and −0.01 • C (R 2 = 0.05; p = 0.50) per year, respectively.The observed LST trends resulted in a total LST increase of 1.05 and 1.56 • C in the morning (10:30) and afternoon (13:30), respectively, and a total decrease of the LST of 0.3 • C (22:30) and 0.12 • C (01:30) at night over the period from 2000 to 2015 in Jambi.
To separate the effect of land use change from global climate warming, we used a site constantly covered by forest over that period (from the forest sites we used in this study) as a reference that was not directly affected by land cover changes.That site showed small changes in LST than the entire province: only the mean morning LST (10:30) had a significant but small trend with an increase of 0.03 • C per year (R 2 = 0.21, p < 0.05) resulting in a total LST increase of 0.45 • C between 2000 and 2015 (Fig. 5a).This LST warming is much smaller than the overall warming at provincial level of 1.05 • C. The LST time series at other times showed no significant trends: the mean afternoon LST (13:30) increased by −0.05 • C per year (R 2 = 0.01, p = 0.31) (Fig. 5b), and the night and evening LST by 0.01 • C per year (Fig. 5c and  d, p = 0.19 and p = 0.60, respectively).
The mean annual NDVI in Jambi decreased by 0.002 per year, resulting in a total NDVI decrease of 0.03 (R 2 = 0.34; p = 0.01; Fig. 5e).The NDVI of the forest showed a small but not significant increase of 0.001 per year (R 2 = 0.04, p = 0.23) (Fig. 5e) fluctuating around an NDVI of 0.84.

Landsat LST compared to MODIS LST
In our study we retrieved the surface temperature from a Landsat image and compared this with MODIS LST.Our results showed a good agreement between both LSTs (Fig. 3), which is comparable to other studies and thus gives confidence in our analysis.Bindhu et al. (2013) found also a close relationship between MODIS LST and Landsat LST by using the same aggregation resampling technique as our method and found a R 2 of 0.90, a slope of 0.90 and an intercept of 25.8 • C for LST, compared to our R 2 of 0.8, slope of 1.35 and intercept of −11.58 • C (Fig. 3).Zhang and He (2013) validated Landsat LST with MODIS LST and also found good agreements (RMSD 0.71-1.87• C) between the two sensors, whereas we found a RMSE of 1.71 • C. Nevertheless, there still are differences and slope versatility between the two satellite sources.These differences are typically caused by differences between the MODIS and Landsat sensors in terms of (a) different sensor properties e.g.spatial and radiometric resolution and sensor calibration; (b) georeferencing and differences in atmospheric corrections (Li et al., 2004); and (c) emissivity corrections; i.e. the use of approximate equations to derive the emissivity from the NDVI from Landsat's red and NIR bands.Li et al. (2004) and Vlassova et al. (2014)   sat, which they attributed to the delay of 15 min in acquisition time between MODIS and Landsat.MODIS LST is measured 15 min later and our results showed that MODIS LSTs were indeed higher than Landsat LST.A comparison of MODIS LST with locally measured canopy surface temperatures during the overpass time of MODIS also showed agreement (Supplement S7, Fig. S7.1).The slope was possibly related to differences in instrumentation and emissivity corrections and to scale issues; this comparison could corroborate the quality check of MODIS LST.
As the MODIS LST product is proven to be accurate within 1 • C (Silvério et al., 2015;Wan et al., 2004) and has been intensively validated, the use of MODIS LST was a proper way to assess the quality of our Landsat LST.
The errors from the different sources (such as atmospheric correction, emissivity correction, resampling Landsat to MODIS resolution) are difficult to quantify.When we tested the impact of atmospheric correction and emissivity errors on the LST from Landsat retrieval we found that (a) the overall patterns across different land use types did not change, (b) emissivity was the most important factor, although the effects on LST retrieval were small, and (c) errors related to atmospheric correction parameters were small because there were minor differences between the default at-mospheric correction (ATCOR) parameters and the ATCOR parameters derived with actual local conditions (RH, air pressure and air temperature).Following the method of Coll et al. (2009) and Jiang et al. (2015) we show that the use of the online atmospheric correction parameter calculator is a good option provided that RH, air temperature and air pressure measurements are available.We additionally compared locally measured air temperatures with MODIS air temperature and found a good agreement (Supplement S8, Fig. S8.1), which served as a verification that we used a correct air temperature for the atmospheric correction parameter calculator.
Overall, our comparison of Landsat LST with MODIS LST and against ground observations suggests that we are able to retrieve meaningful spatial and temporal patterns of LST in the Jambi province.

LST patterns across different land use and land cover (LULC) types
The land cover types in our study covered a range of land surface types that develop after forest conversion.This is the first study in this region that includes oil palm and rubber as land use types that develop after forest conversion.The coolest temperatures were at the vegetated land cover types while the warmest surface temperatures were on the non-vegetated surface types like urban areas and bare land.Interestingly, the oil palm and rubber plantations were only slightly warmer than the forests whereas the young oil palm plantations had clearly higher LST than the other vegetated surfaces.For other parts of the world, Lim et al. (2005Lim et al. ( , 2008)), Fall et al. (2010) and Weng et al. (2004) also observed cooler temperatures for forests and the highest surface temperatures for barren and urban areas.
In Indonesia, land transformation is often not instantaneous from forest to oil palm or rubber plantation but can be associated with several years of bare or abandoned land in-between (Sheil et al., 2009).Oil palm plantations typically have a rotation cycle of 25 years, resulting in repeating patterns with young plantations (Dislich et al., 2016).Given the large LST differences between forests and bare soils or young oil palm plantations that we observed, a substantial warming effect of land transformation at regional scale is expected.

Drivers of local differences between different land cover types
All the land cover types (except acacia plantation forests) had a higher albedo than forest, indicating that these land cover types absorbed less incoming solar radiation than forests.Nevertheless, these land cover types were warmer than forests, suggesting that the albedo was not the dominant variable explaining the LST.Indeed, the statistical analysis showed that ET-LST had a higher correlation than albedo-LST.The ETs were significant, highlighting that despite their higher albedo, all land cover types had higher LSTs than forests related to lower ET rates than forests.In contrast, forests that absorb more solar radiation because of the lower albedo have lower LST because of the higher ET they exhibit, hereby identifying evaporative cooling as the main determinant of regulating the surface temperature of all vegetation cover types (Li et al., 2015).Both observational and modeling studies carried out in other geographic regions and with other trajectories support our observations.Observational studies in the Amazonia by Lawrence and Vandecar (2015) on the conversion of natural vegetation to crop or pasture land showed a surface warming effect.Salazar et al. (2015) provided additional evidence that conversion of forest to other types of land use in the Amazonia caused significant reductions in precipitation and increases in surface temperatures.Alkama and Cescatti (2016) and earlier studies by Loarie et al. (2011a, b) showed that tropical deforestation may increase the LST.Croplands in the Amazonian regions were also warmer than forests through the reduction of ET (Ban-Weiss et al., 2011;Feddema et al., 2005) and that the climatic response strongly depends on changes in energy fluxes rather than on albedo changes (Loarie et al., 2011a, b).A study by Silvério et al. (2015) indeed found that tropical deforestation changes the surface energy balance and water cycle and that the magnitude of the change strongly depends on the land uses that follow deforestation.They found that the LST was 6.4 • C higher over croplands and 4.3 • C higher over pasture lands compared to the forest they replaced, as a consequence of energy balance shifts.Ban-Weiss et al. (2011) and Davin and de Noblet-Ducoudré (2010) added that in addition to the reduction of ET, the reduction of surface roughness most likely enhanced the substantial local warming.
Also for non-Amazonian regions, the replacement of forests by crops caused changes comparable with our observations.In temperate Argentina, Houspanossian et al. (2013) found that the replacement of dry forests by crops resulted in an increase of albedo but still forests exhibited cooler canopies than croplands.The cooler canopies were a result of a higher aerodynamic conductance that enhanced the capacity of tree canopies to dissipate heat into the atmosphere and to both latent and sensible heat fluxes operating simultaneously to cool forest canopies.
In a global analysis Li et al. (2015) showed that tropical forests generally have a low albedo, but still the net energy gain caused by solar energy absorption is offset by a greater latent heat loss via higher ET and that in the tropical forests the high ET cooling completely offsets the albedo warming.For China, this cooling effect was also shown by Peng et al. (2014), who compared LST, albedo and ET of plantation forests, grassland and cropland with forests.
Using NDVI as an indicator of vegetation abundance Weng et al. (2004) (for the US) and Yue et al. (2007) (for China) found that areas with a high mean NDVI had a lower LST than areas with a low mean NDVI, therefore suggesting that vegetation abundance is an important factor in controlling the LST through higher ET rates.Our result support their assumptions by showing the high correlation between NDVI and LST and between ET and LST.
Our findings are also supported by modeling studies.Beltrán-Przekurat et al. (2012) found for the southern Amazon that conversion of wooded vegetation to soybean plantations caused an increase of the LST due to decreased latent heat and increased sensible heat fluxes.Climate models also show the same warming trends and land surface modeling also projects an increase in surface temperatures following deforestation in the Brazilian Cerrado (Beltrán-Przekurat et al., 2012;Loarie et al., 2011b).In a global analysis, Pongratz et al. (2006) showed a LST increase of forest to cropland or pasture transitions, which was driven by a reduced roughness length and an increased aerodynamic resistance, and that the temperature response is intensified in forest to clear/bare land transitions (1.2-1.7 • C increase).Similar to observational studies, the modeling results of Bathiany et al. (2010) show that ET is the main driver of temperature changes in tropical land areas.
To understand the effects of deforestation on biophysical variables in Indonesia, our study identifies the following mechanisms: (a) reduction of ET decreases surface cooling, (b) reduced surface roughness reduces air mixing in the surface layer and thus vertical heat fluxes, (c) changes in albedo change the net radiation and (d) changes in energy partitioning in sensible and latent heat and heat storage.The effect is an increase of the mean temperatures that leads to warming effects in all tropical climatic zones (Alkama and Cescatti, 2016).We point here that our study included a ground heat flux but did not take into account the storage of heat in the soil and the release of stored heat out of the soil during the daily cycle, and the Landsat satellite image was obtained under cloud-free conditions with high shortwave radiation input and low fraction of diffuse radiation.Therefore, the LST retrieved on cloud-free days might be overestimated compared to cloudy days, as the differences in LST between land uses are supposed to be lower when diffuse radiation increases.
Our study is the first to include the oil palm and rubber expansion in Indonesia.In Indonesia, smallholders take 40 % of the land under oil palm cultivation for their account (Dislich et al., 2016).Because the landscape in Jambi is characterized by a small-scale smallholder-dominated mosaic, including rubber and oil palm monocultures (Clough et al., 2016), studies using medium-to coarse-resolution data are not able to capture the small-scale changes and processes at the small-scale level.By using high-resolution Landsat data we were also able to include the effects of land use change on biophysical variables and the underlying processes of the small-scale holder agriculture.

Effects of land use change on the provincial surface temperature in the past decades
The increases of the mean surface temperature in Jambi were stronger during the morning (10:30) and afternoon (13:30) than during the evening (22:30) and night (01:30).Given that our results show a decrease of the NDVI in the same period, this suggests that the observed increased trend of the day time LST can be attributed to the land cover changes that occurred.Our assumption that the observed decreasing NDVI trend is caused by land conversions is supported by two different studies which reported that in Jambi, between 2000 and 2011 (Drescher et al., 2016) and between 2000 and 2013 (Clough et al., 2016), the forest area decreased and that the largest increases were for rubber, oil palm and agricultural and tree crop areas.The class "other land use types", which includes urban areas, showed a minor increase (around 1 %), suggesting that the decrease in NDVI was most likely caused by forest cover loss and not by urban expansion (Table S9).
The same observations on LULC change in Indonesia were also done by Lee et al. (2011), Margono et al. (2012, 2014) and Luskin et al. (2014).Luskin et al. (2014) showed that in Jambi, during the period 2000-2010, forests decreased by 17 % while oil palm and rubber area increased by 85 and 19 %, respectively.Given these trends in LULC changes, the observed LST trends were most likely caused by gradual decrease of forest cover loss at the expense of agriculture and croplands.Our assumptions are supported by findings of Silvério et al. (2015), Costa et al. (2007), Oliveira et al. (2013), Spracklen et al. (2012) and Salazar et al. (2015) that indicate that land use transitions in deforested areas likely have a strong influence on regional climate.Alkama and Cescatti (2016) show that biophysical effects of forest cover changes can substantially affect the local climate by altering the average temperature, which is consistent with our observations and can be related to the observed land use change in the Jambi province.As Indonesia has undergone high rates of forest cover loss from 2000 to 2012 (Margono et al., 2014), these findings support our assumptions that the observed LST increase in the Jambi province was most likely caused by the observed land use changes.
To separate the effect of global warming from land-usechange-induced warming, we considered areas with permanent and large enough forests as a reference where changes are mainly because of global warming.We find that LST of forests shows either no significant trends (at 01:30, 22:30, 01:30) or just a clearly smaller increase of 0.03 • C per year at 10:30.The difference between the LST trend of the province and of the forest at 10:30 was 0.04 • C per year, resulting in a LST of 0.6 • C between the province and forest in the period 2000 and 2015.We point out that our MODIS analysis has a larger proportion of data from the dry season compared from the wet season, as there were more cloud-free conditions during the dry season.Thus, our reported warming effect reflects cloud-free conditions.During cloudy conditions, particularly in the wet season, the warming effect is expected to be lower.A seasonality analysis showed that the relationships in the dry season are stronger than for the wet season (see Supplement S10, Fig. S10.1), which suggests that the warming is more pronounced during the dry season compared to the wet season, which is reasonable as we have more incoming radiation during the dry season.
With the warming effects we found between forest and other land cover types ( LST, Fig. 4a) and the observed land cover changes by Clough et al. (2016) and Drescher et al. (2016) (Supplement S9, Tables S9.1 and S9.2), we estimated the contribution of all land cover types (except forest) to the LST of the province between 2000 and 2015 to be 0.51 • C out of the observed 0.6 • C, which also supports our assumption that the LST increase in Jambi was for 85 % driven by land cover changes (see Supplement S9, Tables S9.1 and S9.2: Land use change analysis), with clearcut areas having a large contribution as they have the largest warming effect.
The observed small but significant increase in LST of forests of 0.03 • C per year at 10:30 reflects a LST change independent of land cover changes, as the forest remained unchanged over that time period.A potential driver of that LST increase is the general global air temperature trend because of changes in radiative forcing or border effects (advection from warmer land uses), which is similar to the 1994-2014 time series analysis of Kayet et al. (2016), who showed a LST increase for all land cover types ranging from wasted land, agriculture land, open forest, dense forest, water bodies and built-up areas.
The observed trends of the provincial air temperature (Fig. 5f) were significant, suggesting that a general warming due to global and regional effects contributes to the observed warming at the provincial level during day and nighttime, but that it is smaller than the land-cover-change-induced effects (Supplement S9, Tables S9.1 and S9.2) at the provincial level (Fig. 5a and b).
In our long-term analysis on the regional effects of land use change we observed an increase in the mean LST and mean air temperature in the 2000-2015 period, concurrent with a decrease of the NDVI.The warming observed from MODIS LST data and from the air temperature obtained from the independent ERA-Interim reanalysis in the Jambi province are most likely caused by the observed decrease of the forest area and an increase of oil palm, rubber and other cash crop areas in the same period, with other effects such as radiative forcing changes and additional natural effects playing a smaller role.Given the plan of the Indonesian government to substantially expand oil palm production with a projected additional demand of 1 to 28 Mha in 2020 (Wicke et al., 2011), the strong warming effect we show for Jambi may serve as an indication of future LST changes for other regions of Indonesia that will undergo land transformations towards oil palm plantations.
A recent study by Tölle et al. (2017) showed that for Southeast Asia, land use change at large scale may not only increase surface temperature but also impact other aspects of local and regional weather and climate, including in regions remote from the original landscape disturbance.Their results also indicate that land clearings can amplify the response to climatic extreme events such as El Niño-Southern Oscillation.The observed effects of land use change on the biophysical variables may have implications for ecosystem services in the Jambi province beyond a pure warming effect.The high precipitation in this region in combination with the reduced vegetation cover of bare land and young oil palm plantations impose risks of soil erosion caused by surface run off.Less water infiltration into the soil, thereby decreasing the soil water storage, may lead to low water availability in the dry season (Dislich et al., 2016;Merten et al., 2016).High surface temperatures in combination with low water availability may make the vegetation and the surroundings more vulnerable to fires.

Conclusions
In summary, we studied the effects of land use and land cover changes on the surface biophysical variables in Jambi and explained the underlying mechanisms of the surface temperature regulation.We derived biophysical variables from satellite data, analyzed the biophysical impacts of deforestation and on a local scale we found a general warming effect after forests are transformed to cash or tree croplands (oil palm, rubber, acacia) in the Jambi province of Sumatra.The warming effect after forest conversion results from the reduced evaporative cooling, which was identified as the main determinant of regulating the surface temperature.On a regional scale, we saw that the effects of land cover changes are reflected back in changes of the LST, NDVI and air temperature in Jambi.The warming effect induced by land cover change clearly exceeded the global warming effect.Understanding the effects of land cover change on the biophysical variables may support policies regarding conservation of the existing forests, planning and expansion of the oil palm plantations and possible afforestation measures.
Data availability.Data are available upon request from the corresponding author.MODIS and Landsat satellite data are distributed by the Land Processes Distributed Active Archive Center (LP DAAC), located at USGS/Earth Resources Observation and Science USGS/EROS Centre, Sioux Falls, SD (http://lpdaac.usgs.gov),and are freely available and accessible.ERA-Interim data are provided by European Centre for Medium-Range Weather Forecasts (ECMWF) and available at http://apps.ecmwf.int/datasets/data/interim-full-daily/levtype=sfc/.
Author contributions.CRS conducted the research, fieldwork and analysis and prepared the manuscript, which was reviewed by GlM, TJ, AM, OR and AK.AM and AK provided the meteorological data.
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.Geographic location of the study area.Jambi province on the Sumatran Island of Indonesia (a, b).The background of the map (c) is a digital elevation model, showing that the plots are located in the lowlands of the Jambi province.The large rectangles are the 28 different land cover types; the small squares are the locations of the 49 small plots of the 7 different land cover types.The land cover types are abbreviated as CLC (clear-cut land), UB (urban area), YOP (young oil palm plantation), MOP (mature oil palm plantation), PF (acacia plantation forest), RU (rubber plantation) and FO (forest).
Landsat 7 ETM+ VIS/TIR 30 m resolution surface reflectance image with low cloud cover, acquired at 10:13 (local time) on 19 June 2013 covering the lowland area of the Jambi province (path 125, row 61), was used in this study.Like all Landsat 7 ETM+ images acquired after 31 May 2003, the image we used was affected by a scan line error causing a data loss of about 22 % (http:// landsat.usgs.gov/products_slcoffbackground.php).Most selected plots were located in the center of the image and thus not affected by the data loss; e.g. the forest plots located at the edges of the scan line error zone faced minimal data loss because they were large enough.We also downloaded the tile h28v09 of the MODIS Terra (MOD) and Aqua (MYD) daily 1 km Land Surface Temperature and Emissivity products (MOD11A1 and MYD11A1 Collection 5) and MODIS 16-day 500 m vegetation indices NDVI/EVI product (MOD13A1 Collection 5) from 5 March 2000 to 31 December 2015 for Terra data and from 8 July 2002 to 31 December 2015 for Aqua data.We downloaded other supporting satellite data such as the MODIS Atmospheric Profile product (MOD07_L2) and the MODIS Geolocation product (MOD03).All MODIS data were reprojected to WGS84, UTM zone 48 south with the MODIS Reprojection Tool (MRT).The quality of the MODIS data was examined with the provided quality flags and only pixels with the highest-quality flag were used in the analysis.2.4 Retrieval of biophysical variables from Landsat 7 ETM+ VIS/TIR images NDVI NDVI was derived from the reflectances corrected for atmospheric effects in the red (ρRED, band 3 Landsat 7 ETM+) and near-infrared (ρNIR, band 4 Landsat 7 ETM+) bands, with NDVI = ρNIR − ρRED ρNIR + ρRED .

Figure 2 .
Figure 2. MODIS LST image (a) compared with Landsat LST image (b).Cloud cover and cloud shadow cover resulted in data gaps (no data).The difference in acquisition time between the images is 15 min.The square in the MODIS image is the area that is covered by the Landsat tile (path 125, row 61).Both satellite images were acquired on 19 June 2013.

Figure 3 .
Figure 3. Average surface temperature (LST) and standard deviation (SD) of seven land cover types derived from a Landsat thermal image compared with the mean and SD of MODIS LST.The dashed line is the theoretical 1 : 1 line, and the solid lines are the linear model type 2 regression line (black) and the confidence limits of the regression line (red).The Landsat and MODIS images were acquired on 19 June 2013, at 10:13 and 10:30 local time, respectively.Landsat pixels (30 m) were resampled to MODIS pixel resolution (926 m) to make a pixel-to-pixel comparison between the two sources possible.

Figure 5 .
Figure 5. Mean annual LST (a-d), mean annual NDVI (e) and mean annual air temperature trends (f) in the Jambi province between 2000 and 2015 derived from MODIS LST (a 10:30, b 13:30, c 22:30 and d 01:30, local time), MODIS NDVI and ERA-Interim daily air temperature (01:00 and 13:00, local time) data sets, respectively.Grey-shaded areas are the confidence intervals of the means; blue-shaded areas are the confidence intervals of the regression lines.MODIS LST time series for 13:30 and 01:30 were available from the middle of 2002; for this reason we used the complete years from 2003 to 2015.

Table 1 .
Statistical analysis between biophysical variables (albedo (α), NDVI and ET) and spectral radiance band (L6), corrected thermal band (R c ) and Landsat land surface temperature (LST).* * p = 2 × 10 −16 .LM is multiple linear regression analysis between LST (or L6 or R c ) and three biophysical variables: albedo (α), NDVI and ET.ρ is the correlation coefficient.R 2 is the R square of the components.β is the regression coefficient of the component.Stand.β is the standardized β.Model fit (R 2 ) is the overall model fit of the multiple linear regression. *