Spatiotemporal Distribution of Human–Elephant Conflict in Eastern Thailand: A Model-Based Assessment using News Reports and Remotely Sensed Data

In Thailand, crop depredation by wild elephants intensified, impacting the quality of life of local communities and long-term conservation of wild elephant populations. Yet, fewer studies explore the landscape-scale spatiotemporal distribution of human–elephant conflict (HEC). In this study, we modeled the potential HEC distribution in ten provinces adjacent to protected areas in Eastern Thailand from 2009 to 2018. We applied the time-calibrated maximum entropy method and modeled the relative probability of HEC in varying scenarios of resource suitability and direct human pressure in wet and dry seasons. The environmental dynamic over the 10-year period was represented by remotely sensed vegetation, meteorological drought, topographical, and human-pressure data. Results were categorized in HEC zones using the proposed two-dimensional conflict matrix. Logistic regression was applied to determine the relevant contribution of each scenario. The results showed that although HEC probability varied across seasons, overall HEC-prone areas expanded in all provinces from 2009 to 2018. The largest HEC areas were estimated during dry seasons with Chantaburi, Chonburi, Nakhon Ratchasima, and Rayong provinces being the HEC hotspots.However, the HEC potential was reduced during severe and prolonged droughts caused by El Nino events. Direct human pressure caused a more gradual increase of HEC probability around protected areas. On the other hand, resource suitability showed large variation across seasons. We recommend zone-dependent management actions towards a fine-balance between human development and the conservation of wild elephants.


Introduction
Although wild Asian elephants (Elephas maximus) are listed as endangered species under the IUCN (International Union for the Conservation of Nature) Red List of Threatened Species, their conservation is hampered by humans as a reaction to crop depredation caused by elephants [1]. In Asia, approx. 10%-15% of total agricultural output can be damaged by wild elephants, which threaten human security and well-being [2]. India hosts the largest population of wild Asian elephants: each year, 400 people and 100 elephants die as a result of human-elephant conflict (HEC) [3]. Understanding the HEC phenomenon is critical to both conservation success and the livelihood of human communities in close proximity to wild elephant habitats.
Thailand is estimated to have 3000-3500 wild elephants in 68 areas, 41 of which are facing HEC, commonly in the form of crop depredation [4]. Historically, elephants have been recorded inside The aim of this study was to bridge the knowledge gap on the physical factors that potentially govern the spatial distribution of HEC in Eastern Thailand. Given the large landscape extent and the dynamic spatiotemporal variation of environmental and physical factors, we utilized remotely sensed satellite data and quantified the spatiotemporal HEC distribution over a 10-year period. Our specific objectives were to (i) model the potential spatial distribution of seasonal HEC from 2009 to 2018 with the use of time-calibrated SDMs, (ii) identify and distinguish the contribution over time of important modeling factors (resource suitability and direct human disturbance), and (iii) prioritize the areas that require targeted management and increased intervention.

Study Area
Our study was carried out in two forest-dominated areas of Eastern Thailand (Figure 1) covering eight eastern and two north-eastern provinces. The area has a tropical monsoon climate. The monsoon season occurs from mid-May to mid-October with an average rainfall of 1400 mm, while the area during dry season receives about 400 mm of rainfall [38,39]. The region has nine national parks (NP) and wildlife sanctuaries (WS) hosting elephants [5]. Khao Angruenai-WS, for example, experiences high density of wild elephants, approx. 0.2 elephant/km 2 [40]. In addition, a constantly low elevation in the central region enabled elephants to easily disperse into agricultural land. Consequently, this area is suspected to be a HEC hotspot. Agriculture is the dominant land cover. The five most important crops in planting areas are rice, cassava, rubber plant, sugarcane, and maize. Orchards and plantations commonly spread out in the southern areas. To limit the modeling boundary to only those potentially accessible by elephants, a 20-km buffer was created from the boundary of protected areas and village location with reported HEC. The seasonal models were set according to the monsoon pattern, with May to October as the wet season and November to April of the following year as the dry season.

Data
HEC occurrences were collected from online news reported between 2014 to 2018. The environmental predictors were prepared under the same time period. Based on this data, we modeled the spatial distribution of HEC across the study area for the period 2009-2018 and developed maps of HEC category. The flow chart of this study is shown in Figure 2.

HEC Occurrence Data
Until March 2019 when elephant-induced damages were first included in farmers' insurance schemes, HEC was neither compensated, nor insured by the Thai government [11]. Consequently, official records were not consistently maintained across protected areas. Although elephants' locations outside of protected areas were sometimes documented, the presence of elephants is not always equivalent to HEC. Reporting from news sources usually happens when negative outcomes occur, and this better reflects HEC occurrences. Therefore, HEC incidences were retrieved from online news sources. 'Wild elephants' in Thai language was used as a search keyword from the News section of Google Search Engine. A customized time period was set between 2014 and 2018. Each search output was investigated manually to exclude duplicated reports of the same incident.
The news reports, however, did not mention the precise locations of HEC occurrence but only the village names. To overcome this lack of exact occurrence locations, we simulated the occurrence locations using a conditional random sampling method. The sampling boundaries were restricted within a 3-km buffer around the center of each mentioned village, excluding areas that fall within the protected areas, large water bodies (e.g., reservoirs), and major road networks. We excluded locations with the aforementioned features because HECs are unlikely to occur within them. The numbers of random occurrence points were generated according to the numbers of damage incidents reported within each village. A total of 124 incidents occurred in the wet season; a combination of 7, 12, 14, 31, and 60 incidents from 2014 to 2018 respectively. The dry season had 122 reports in total; 5, 20, 20, 20, and 57 of which occurred respectively in the same period of time.
Five sets of random occurrences were generated. A Wilcoxon-Mann-Whitney test [41] was performed to compare the distribution profile of each independent variable to that of the other four sets of simulated occurrence records. A p-value of over 0.05 indicated no significant difference between each set. We then selected one set for our model constructions.

Predictor Variables
Thirteen variables were analyzed. We grouped the predictors into two scenarios: (i) resource suitability and (ii) direct human pressure.
Resource suitability comprised of vegetation productivity (Enhanced Vegetation Index-EVI), seasonal vegetation changes (EVI slope, and EVI standard deviation), landscape composition (EVI homogeneity), meteorological drought condition (Keeetch-Byram Drought Index), refuge locations (Forest percent cover, Distance to forest), and topographic condition (Terrain Roughness Index). Direct human pressure included distance to lit-up area, to main roads, to protected habitats, and human population density. Indirect human pressures, such as sociopolitical factors, were not considered in this study. All the predictors were re-projected to the WGS 84/UTM zone 47N (EPSG:32647) and resampled to a 500-m resolution using bilinear interpolation. Pre-processing was performed using Google Earth Engine [42] and R version 3.5.3 [43]. Table 1 shows each variable together with the data source, the original resolution, and the temporal period used. Normalized Difference Vegetation Index (NDVI) was found to be an effective proxy of forage availability [44]. Dispersal of African elephants was shown to coincide with the greening-up measured by NDVI [45]. For this study, we utilized the Enhanced Vegetation Index (EVI). Despite being similar to NDVI, EVI improved saturation in high biomass regions, corrected for aerosol influence, and reduced noise from soil background [46]. Following [47], EVI was calculated from MODerate Resolution Imaging Spectroradiometer (MODIS) Terra product (MOD09A1) as: where ρ NIR , ρ Red , and ρ Blue represent the reluctance of the near-infrared, red, and blue bands respectively. Only the pixels under clear cloud state and no cloud shadow were used. We first calculated the monthly median of EVI for each month from 2009 to 2018. The missing monthly pixels were filled using the 10-year averaged EVI value in the same pixel location of the same month. From the monthly EVI data, we calculated mean EVI for each season which represents vegetation productivity. Next, the EVI slope variable, representing the rate of change in vegetation condition (e.g., crop senescence), was calculated by applying pixel-wise linear regressions over the monthly EVI within each season.
A standard deviation of the monthly EVI values within each season was calculated next which represents fluctuation in vegetation dynamic. Lastly, a spatial homogeneity of EVI was generated using the Gray Level Co-Occurrence Matrix [48]. These EVI variables can also be linked to different characteristics of land cover types. For example, a high EVI and a EVI slope near zero usually associate with tropical forest land cover ( Figure S1). Drought influences surface water availability and vegetation quality, which govern elephant habitat use [28]. The Keetch-Byram Drought Index (KBDI) estimates the dryness of soil layers. The KBDI product was computed using the precipitation data derived from the Global Satellite Mapping of Precipitation (GSMap) and land surface temperature (LST) data from Multi-functional Transport Satellite (MTSAT) [49]. The value of KBDI ranges from 0 (no moisture deficit) to 800 (extreme drought). The daily data from 2009 to 2018 was averaged by season. Additionally, wild elephants were observed to move toward inland areas during the dry season as waterhole in coastal regions dried up [29]. To capture accessibility to water, locations of surface water were obtained from the monthly historical Landsat Global Water Surface Product [50]. Within a single year, pixels detected with water for at least 3 months were marked as water and Euclidean distance to them were calculated.
Forest is considered a natural habitat and represents a potential refuge location. Forest land cover classes from the MODIS land-cover product (MCD12Q1) was used. Since the study area was dominated by dry evergreen forest (90% of all forest classes), we reclassified all forest types to a single land cover class. According to an interview with park rangers conducted by the authors, 6 km was suggested to be a one-way distance traveled by wild elephants between patches of the forest outside the protected areas. Two variable were calculated, a mean Euclidean distance from each pixel to forest and a percentage of forest cover within 6-km buffer around each pixel. Lastly, we calculated terrain ruggedness index (TRI) from the Shuttle Radar Topography Mission data (SRTM) [51]. The TRI represents the relative change in elevation from a center cell and eight surrounding cells. A higher TRI value indicates more rugged areas.
Direct human disturbance was measured based on human population density, as well as Euclidean distance to protected habitats, to main roads, and to lit-up areas. The pixels detected with light or the lit-up pixels were computed from satellite-derived night time light data. The Defense Meteorological Satellite Program's Operational Linescan System (OLS) and the Suomi National Polar-orbiting Partnership satellite's Visible Infrared Imaging Radiometer Suite (VIIRS) were the main sources of night-time light product. Calibration among OLS sensors, as well as between OLS and VIIRS, was necessary. We applied a second-order regression model from [52] for the calibration among OLS sensors. For OLS and VIIRS inter-calibration, we first created VIIRS annual composite following [53] and then applied a combination of power function and Gaussian low pass filter [54]. The pixels with digital numbers of over 20 were used to calculate Euclidean distance. The density of the human population also influences alteration of landscape and intensity of anthropocentric activities, which may not be captured by the night-time lights. Hence, mean human population density was computed using yearly estimations from LandScan.

Bias Correction
HEC incidences from online news sources are opportunistically collected and not randomly sampled. Such datasets often contain sampling bias wherein more reporting are made from easily accessible locations or well-known hotspots [55]. With sampling bias, it is hard to determine whether occurrences were reported due to preferable conditions in that locations or concentration of search effort. When relative search effort across the landscape is known, sampling bias can be directly modeled and provided as prior distribution during SDM construction [56,57]. Alternatively, the effect from sampling bias can be partially accounted for by subsampling the training dataset or adjusting the background selection [55,58,59]. Due to low occurrences in our study, we applied background selection method which nullifies bias by generating a similar bias in the background [57]. Bias grids were produced by deriving a Gaussian kernel density map of the village locations weighted by the average number of duplicated reports within each village. The bias values were re-scaled from 1 to 20, following [58] to avoid extreme values, and used as probability in sampling background points. We sampled a total of 10,000 points, a combination of 2000 each year from 2014 to 2018. The generated background points were later used as pseudo-absences in model construction.

Maximum Entropy Modeling
The Maximum Entropy algorithm from MaxEnt [60] was used. MaxEnt is a machine-learning technique that estimates the unknown distribution of suitability by contrasting the values of predictors at occurrence locations with the overall distribution of these predictors [56]. A detailed explanation and related equations can be found in [60]. MaxEnt had shown a high performance even with few occurrence records and was least affected by errors of occurrence location [56]. It also outperformed other methods [61]. All our models were constructed using dismo package in R with MaxEnt 3.3.4 version [62]. The logistic link function was used to derive a relative probability of potential HEC occurrence ranging between zero (low probability) and one (high probability) [60].
A time-calibrated method [63] was applied in which each occurrence point was matched with environmental predictors from the relevant season during which HEC was reported. This resulted in time-independent models which allowed comparability across the study period. MaxEnt requires background points as pseudo-absent. The 10,000 background samples previously generated in Section 2.3.1 were used.
Prior to model construction, multicollinearity among predictors was evaluated. Variables that had Variance Inflation Factors (VIF) greater than 10 and high Pearson correlation (−0.75 < r < 0.75) were removed. Since feature classes and regularization multiplier (RM) impacted modeling results [56], parameter optimization was conducted using EMNeval package in R [64]. Product and Threshold features were excluded in our models. Product-feature tends to over-fit and complicates interpretation of variable responses [65]. Threshold-feature should be used when a drastic cut-off exists in species' response to environmental factors, but no such cut-off has been identified for Asian elephants. Therefore, only Linear/Quadratic/Hinge combinations were selected and k-fold cross-validation was performed with RM value from 0.5 to 5 at 0.5 increments. Akaike Information Criterion (AIC) was used for optimal parameters selection. Other settings were left with default values which included 500 iteration maximum and convergence thresholds.
After optimal parameters were identified, the models were constructed using k-fold cross validation (k = 5) and evaluated with Receiver Operating Characters (ROC) with average Area Under the Curve (AUC) from all replicas. In addition, a jackknife test was used to identify important predictors. Responses for each variable were also generated. A total of four models were constructed, one model for each season under the two high-level scenarios (resource suitability and direct human pressure). We identified the differences between environmental predictors from each season to those used for model construction using Multivariate Environmental Similarity Surface (MESS) and limiting factors [58]. The negative MESS score indicated a novel condition in variables used for prediction which implies possible uncertainty. We then estimated relative probability of HEC across the landscape for 20 seasons by applying our constructed models on the predictors from 2009-2018.

Conflict Classification and Analysis
The probability of HEC occurrence for each high-level group was then categorized in three classes (High, Low, Very Low). Two thresholds were used, (i) 10th percentile of presence locations and (ii) maximum training sensitivity plus specificity (maxSS). The first threshold allowed omission of 10 percent of occurrences which reduces sensitivity to extreme localities [66], while the second threshold was evaluated as an effective threshold value for presence-only modeling [67]. Probability lower than the first threshold was set as Very Low class. We then applied the second threshold where the probability lower than maxSS was set as Low, and those higher or equal to maxSS was set as High.
Interpretation of resource suitability is straight forward in which high HEC probability occurred in a more suitable condition. Conversely, the high HEC probability of human pressure captured the disturbance level in which conflict peaked. In reality, high human disturbance beyond the peak level existed, but likely restricted occurrence of elephants resulting in low predicted HEC probability.
Each classified maps from different scenario in the season were then overlaid into a two-dimensional HEC categorical map (Figure 2). These categorical classification contained Avoid matrix (at least one very low class from either scenario), Rare conflict (low resource suitability and low human pressure), Low conflict (high resource suitability but low human pressure), Likely conflict (low resource suitability and high human pressure), and High conflict (high resource suitability with high human pressure). By using two-dimensional classification, we can identify two main key management-relevant actions related to each group of factors. First, management actions associated with resource suitability are linked to natural resource and land management (e.g., land-use policies, establishing elephants corridors) (e.g., [68,69]). Second, HEC occurrences are also governed by level of human disturbances which can be associated with different management actions directed more toward human co-adaptation (e.g., insurance schemes, behavioral adjustment in crop husbandry) (e.g., [70,71]).
We generated two HEC maps (a wet seaon map and a dry season map) for each year during 2009-2018. Areas of different HEC levels were calculated by summing the number of pixels within each category and multiplying that by the pixel size. We calculated affected areas for each map by season from 2009 to 2018 both for the whole region and separately by provinces. The distributions of conflict hotspots, which are the areas repeatedly predicted with the same conflict category across the years, were identified. The change in probability of HEC occurrence under resource suitability and direct human pressure scenario over 10 years was calculated by fitting pixel-wise linear regression on predicted probability from 2009 to 2018. The slope coefficient of the fitted regression indicated the rate and direction of change in HEC probability. The intercept represented the baseline probability in 2009.

Model Performance and Variable Responses
The p-value of all simulated HEC occurrences were greater than 0.19. The mean cross-validated AUCs were 0.81 for resource suitability/wet season, 0.73 for resource suitability/dry season, 0.78 for direct human pressure/wet season, and 0.77 for direct human pressure/wet season. Jackknife analysis for the resource suitability scenario (both seasons) indicated that Forest Percent Cover had the highest predictive contribution. Together with KBDI, EVI slope, and Distance to Forest, these four predictors accounted for over 80% of the models' predictive power. For the dry season model, KBDI had slightly less contribution, while distance to forest edge became more important. Under human pressure scenario, Distance to Protected Habitats was the most influential variable for both season. The next important predictor for the wet season was Human Density, while Distance to Main Roads was for the dry season. Figure 3 shows how the HEC probability changes as each environmental predictor is varied, while keeping all other environmental variables at their average sample value. The response of Forest Percent Cover was similar for both seasons. In the highest and lowest forest densities, lower HEC probability was expected, while higher HEC probability was found in moderate forest densities. For KBDI, higher probability of HEC in wet season occurred at low KBDI, with a continuous reduction after KBDI of around 50. In the dry season, however, HEC probability peaked at intermediate KBDI around 300-400 and decreased slowly. The response of EVI Slope for the wet season indicated that HEC was more likely to occurred where vegetation conditions were changing (diverted from zero), the highest HEC probability occurred when EVI was reducing over the season. In the dry season, however, probability of HEC was higher when EVI was relatively stable (EVI of zero) or increasing slightly (EVI of 0.05), indicating green-up of vegetation. The patterns of EVI slope from both seasons corresponded to the characteristics of EVI slope from forest and savanna land cover ( Figure S1). For human pressure, response of Distance to Protected Habitats, Distance to Main Roads, and Distance to Lit-up Areas had similar characteristics for both seasons. In the dry season, a slower reduction in HEC probability was observed for Distance to Protected Habitat and Distance to Lit-up Areas as distance increased. For Human Density, HEC probability in the wet season reduced as density approached 1000 person/km 2 , but did not affect HEC probability in the dry season. Possible higher tolerance to high human pressure of elephants was captured in dry season.

Distribution of Conflict and Conflict Hotspot
The potential for HEC occurrence was higher during the dry season: High and Low conflict areas were larger and more frequent. In contrast, during the wet season, the Likely and Rare conflict categories were more frequent, suggesting lower HEC potential. The hotspots of High conflict category were concentrated around the south and south-west of Ang Ruenai-WS in Chonburi, Rayong, and Chantaburi provinces (Figure 4). In the north, smaller clusters, especially near the protected areas, was predicted in Nakhon Ratchasima, Nakhon Nayok, and Prachinburi provinces. The high HEC zones shrunk closer to the protected areas in the wet season and mainly located around Khao Chamao Khao Wong-NP at the border between Rayong and Chantaburi, east of Khao Soi Dao-WS in Chantaburi, and northwest of Khao Yai-NP in Nakhon Ratchasima. Additionally, our models estimated that many areas under High conflict category in the dry season changed to Likely conflict category in the wet season. This result implied that such locations have potentially experienced year-round HEC in different levels (e.g., intensity or frequency). Low conflict class was predicted in large areas in the dry season, affecting all provinces. Although less in frequency and intensity, in the wet season coldspots of Low and Rare conflict categories were concentrated around the main roads, with some areas located far from protected areas. During 2009-2018, overall areas of potential conflict were estimated to be increasing as shown in Figure 5. The increasing trend of High HEC was captured in both dry and wet season, although the peak values were lower in the wet season. Potential HEC areas expanded more than double between 2016 and 2017, of which areas with High conflict increased from 2235 to 4306 km 2 and 115 to 2467 km 2 in the dry and wet season respectively. The dry season was dominated by two conflict categories: High and Low. On the other hand, similar trends were presented across all conflict categories in the wet season despite variations of affected areas among the years. For Likely conflict, the wet seasons had a similar increasing trend to that of the High conflict category. However, the dry seasons had relatively stable areas of the Likely and Rare HEC.
Chantaburi was estimated to have the largest areas of HEC, followed by Nakhon Ratchasima (Figure 6). In Chantaburi, large areas of High HEC (~900 km 2 ) were estimated in the dry season from the beginning of the study period. The province showed an increase in overall areas of conflicts, as well as the largest area expansion of High HEC captured in the wet season, from 170 km 2 in 2009 to 689 km 2 in 2018. Nakhon Ratchasima also had large HEC-prone areas, but the High conflict category showed a large increase only from 2014 onward. Similar to Nakhon Ratchasima, Buri-Ram and Chachoengsao were predicted with High HEC from 2014. Except Nakhon Nayok and Trat, all provinces were predicted to have a larger area of High conflict category during the dry season. HEC areas were increased more than double from 2016 to 2017 in Buri-Ram, Chachoengsao, Chantaburi, Nakhon Ratchasima, Prachinburi, Rayong, and Sa Kaeo. On the other hand, a decrease in the areas of HEC was identified in 2010 and 2014-2016 for most provinces.

Drivers of Changes in HEC Probability Over Time
We identified the contribution to changes of HEC by evaluating HEC probability from resource suitability and direct human pressure across the study period. From Figure 7a, HEC probability from direct human pressure scenario generally showed a gradual increasing trend with an exception of a drastic area expansion in both High (2203 to 6503 km 2 ) and Low (4773 to 8983 km 2 ) classes in 2014. This sudden increase was likely caused by lit-up areas increased as a result of improved night-time light sensor started from 2014 onward. For resource suitability scenario, a clear pattern cannot be observed. Hence, variation of predicted HEC category seen among different years were likely due to the dynamic changes in suitable resources. Areas of High and Low probability under resource suitability were reduced over half in 2010 and seemed to continuously decrease from 2012 to 2016. This reduction in HEC areas coincided with the high anomaly of KBDI period in Thailand. Figure 8   In resource suitability maps, areas with orange color corresponded to a high baseline with moderate negative rate of change in HEC probability. This decreasing trend (−0.07 to −0.04 per year) in both wet and dry season was mainly predicted around the edge of the forests. Base on MODIS land cover (Figure 9), the reduction of HEC probability near the forest was due mainly to forest cover increased over the years. According to the predictors' responses in Figure 3, areas with high forest densities and nearer to forest were estimated with lower HEC probability. Positive trends (0.2 to 0.4 per year) was sparsely predicted in both wet and dry seasons on the west-side of Ang Ruenai-WS in Chachoengsao and Chantaburi provinces. Although we cannot specify a reason behind this increase, we observed from MODIS land cover that there was an expansion of the savannas land cover during 2017-2018, as well as increased of forest in those areas (Figure 9b). The high HEC probability of EVI slope corresponded to the characteristics of forest and savanna which may have heightened the predicted probability.  For direct human pressure scenario, large areas of increasing trend occurred in previously low and moderate HEC probability baseline in wet (a pure bright blue) and dry (a dark greenish-blue) season respectively. These areas were located around Ang Ruenai-WS, north-east of Khao Yai-NP, and north of Thablan-NP. Since the variables used under direct human pressure scenario only contain static and annual characteristic (Table 1), the differences observed between seasons were not the result of physical differences. The dissimilarities were due to seasonal differences of variable response that governed HEC prediction. In addition, the southern areas of Ang Ruenai-WS were predicted with constantly high HEC probability (around 0.7) from human pressure in the dry season. The increasing trend within the same areas in the wet season (a rate of 0.07 to 0.10 per year) likely caused a year-round HEC. The large positive trend from direct human pressure, when happen in areas with already high HEC probability predicted under resource suitability, may escalate HEC to a higher category. Lastly, spares areas in orange indicated a decreasing trend in HEC probability. These areas scattered in the west near the main roads. This is due to human population growth. The same was not predicted for the dry season because of the differences in variable response between the two models.

Discussion
Using two-dimensional classification together with time-calibrated SDM on remotely sensed satellite data, we predicted and compared potential HEC distribution in Eastern Thailand across 20 seasons. Overall, our models predicted the occurrence of large HEC-conflict areas during the dry season which decreased both in term of spatial extent and intensity during the wet season. Chantaburi and Nakhon Ratchasima were predicted to have the largest HEC-prone areas. Drought-induced decrease in the distribution of suitable resources (base on KBDI index), resulted to a relevant decrease in the spatial extent of HEC. The high KBDI detected in 2010 and 2014-2016 ( Figure 8) coincided with the El Nino phenomena that caused severe and prolonged drought in Thailand [72]. This caused a decrease in the spatial HEC extent in some provinces, but HEC extent increased in other provinces. Previous studies in arid savannas showed that extreme drought can alter the distribution and abundance of elephants population, leading to mass starvation [73,74]. However, such extreme events is usually not considered during modeling. Considering that dry periods, and their associated extreme drought events, may occur more frequently due to climate change, HEC distribution may become unpredictable, causing critical management implications. Additional field investigation is required to examine whether or not a decrease in the spatial HEC extent in some areas may result to concentrated increase of elephant-induced damages in other locations.
The peak of HEC probability for forest percent cover, a variable with the highest predictive power, was identified at 25%-45%. Forest land cover mainly entailed protected areas. Hence, HEC occurred in areas close to elephant natural habitats. Although high HEC occurrence near protected areas is expected, the peak probability implies that conflict incidents do not always locate directly adjacent to protected parks. Available patches of forest outside of protected areas, such as community forest, may assist in wild elephant dispersal as long as the composition with other land cover provide 25%-45% cover within 6 km. Further field study is necessary to identify the size of forest patches required by elephants outside of protected areas.
HEC hotspots along the southern and western of Ang Ruenai-WS were dominated by savannas land cover, a mixed tree and grass system. The peak in HEC probability as seen from EVI slope coincided with characteristics of MODIS savannas and forest land cover class. By comparing MODIS land cover with land use map from Land Development Department, Ministry of Agriculture and Cooperatives of Thailand, savannas land cover was generally rubber plantations and orchards. Available tree canopies in these land use, despite being sparse, may provide cover for elephants and assist in their movement. Studies in India and China identified proximity to forest edge and high-statue vegetation (e.g., eucalyptus and acacia) as important factors determining elephant occurrences outside of protected areas [25,75,76]. With large continuous extant of savannas within these hotspot, together with high HEC predicted across both season, elephants may already be frequent and even residing permanently in the areas.
Although predictor responses were generally similar to studies from other Asian countries, some variable contributions were different. Distance to water has been identified as an important factor determining elephants distribution in China [25], Indonesia [77], India [78], and Thailand [79]. In this study, however, it was not a prominent predictor. We expected the reason to be the coarse spatial resolution as we re-sampled the data from 30m to 500m. Consequently, small water bodies located in individual farmers' lands might not be captured. Although vegetation index was a good proxy for forage quality in [44], the EVI variable had relatively lower predictive importance in our study. Usefulness of vegetation index depends highly on the types of habitat and the season [80]. In tropical forest, elephant forage abundance was unable to be mapped directly using average value of NDVI [81]. Nevertheless, EVI slope had the highest predictive power among all EVI variables, implying possible importance of vegetation phenology. A study of crop raiding behavior in African elephants identified crop availability and ripening timing as important indicators for predicting crop damage [82]. Therefore, variables related to crop types along with its phenology, which can be detected from remotely sensed satellite data, should be further studied.
Overall our calibrated models had a good to very good predictive power with AUC ranging from 0.73-0.81. Nevertheless, this study still contains limitations and uncertainties. First, although bias correction was applied, we still cannot account for unreported locations where HEC occurred but not reported in the news. Second, different in variable responses were identified between wet and dry season, but the significance between predicted wet and dry HEC distributions remain to be evaluated. Having two separate HEC maps can support effective operational planning (e.g., seasonal patrol routes), but can also cause difficulty for policy-level planning. Future study can evaluate models from key season similar to [35] in which winter season was chosen. Third, additional variables can be included to provide better prediction of HEC occurrence. Besides potential use of cropping pattern and phenology, human tolerance and perception of risk should also be considered. These factors represent the possibility of coexistence between human and wildlife [16]. Fourth, current mitigation efforts have not been included, but are essential as they can alter elephants' access to resources. Previous studies have shown that implementation of physical barriers shift HEC to new locations [83]. Such data on existing mitigation can be incorporated after to identify movement routes and potential corridors. Lastly, future study can include assessment of habitat quality within protected areas. Our models predicted an increase in potential HEC areas after prolong drought during 2014-2016. We cannot infer a conclusion based on our current study. However, extreme events can cause a change to land use and vegetation in the following years, impacting elephant's natural habitat and adjacent agricultural lands. Such information can elucidate root cause of conflict and enhance management decision.
Making informed decision on where to allocate limited resources is crucial for government and conservation organizations alike [84]. A two-dimensional classification approach has been used to identified management-relevant actions in habitat modeling [35][36][37]. Utilizing similar method, we recommended prioritization of HEC-zone dependent management actions. Two groups of management actions are considered, (i) natural resource/land management and (ii) promotion of human adaptation. High HEC zones must receive first priority with parallel emphasis on both land-use policies and human adaptation. Together with HEC-relevant land management, behavioral adaptation of those who live in the areas are important to reduce risky behaviors. In Likely HEC zones, certain land management actions are not necessary (e.g., permanent electric fences), but more focus should be put on community development. For areas with Low HEC category, the extensive change in human behavior may not be needed (e.g crop husbandry), but general knowledge on appropriate actions when encountering with wild elephants are potentially useful. In such areas, land use planning is more important in preventing further escalation of conflict. Lastly, Rare HEC zones were predicted sparsely and far from protected areas (Figure 4). Since these areas were concentrated closer to the main roads, management may focus on the risk of vehicle-elephant-collision. Further field investigation and data collection are necessary to pinpoint appropriate management actions.

Conclusions
This study utilized publicly available dataset and applied time-calibrated SDM with two-dimensional conflict classification to estimate time-series distribution of potential HEC in eastern Thailand. We illustrated the impact of inter-annual and seasonal variation on the distribution of potential HEC, as well as distinguished the contribution between resource suitability and direct human pressure. Overall increasing trend of HEC was predicted from 2009 to 2018, with larger extend in the dry season. Large reduction of conflict areas in 2010 and 2014-2016 were likely explained by server drought due to El Nino events which was captured by KBDI. Our results also suggested that variation in probability and distribution of HEC was due to changes in resource suitability, while a more continuously gradual increase was observed from direct human pressure. Besides identifying HEC response to environmental characteristics, our findings can also support prioritization of conservation resources. We recommended HEC-zone dependent management focus. Parallel emphasis on extensive land management and human co-adaptation should be performed in High HEC-zones. In Likely HEC-zones, more attention should be given to raise safe behaviors for communities. Land use policies in Likely HEC-zones should be strengthened with general awareness of appropriate actions when encountering wild elephants. Rare HEC-zones were scattered close to main roads, hence we recommended investigation to prevent vehicle-elephant collision. Within each zone, more priority can be given to hotspots that estimated with repeated HEC. Lastly, this study highlighted the advantages of satellite-derived variables with high temporal resolution which can capture annual and seasonal variation.
Author Contributions: N.K. and W.T. conceived and designed the study; N.K. performed data curation; N.K. analyzed the data; W.T. supervised and provided resources; N.K wrote the paper; W.T. reviewed and provided revision suggestions. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.