December 2020 - IJG

GIS-Based Landslide Susceptibility Map Verification by its Geotechnical and Geological Characteristics Herianto, P.,* Kunsuwan, B., Mairaing, W., Chalermpornchai, T. and Srisook, W., Division of Geotechnical Engineering, Department of Civil Engineering, Faculty of Engineering at Kamphaeng Saen, Kasetsart University, Kamphaeng Saen Campus, Thailand, E-mail: [email protected] *Corresponding author Abstract The famous 2006 Laplae rainfall-triggered landslide has caused hundreds of million Baht lost. The objective of this paper is to evaluate and to compare present geotechnical and geological conditions between the landslide area and landslide-free area to verify a landslide susceptibility map. Detailed field surveys acquired 14 samples: undisturbed, disturbed soils as well as rocks. Sensitivity and back-analysis were also performed. 35.44% of landslides occurred in meta-shale - greywacke-interbedding Formation, and another 64.56% happened in meta-greywacke - basaltic andesite Formation. The dominant slope of the landslide is in 300-400 slope, with aspects toward South direction, and in a range of <1,000 meters from geological structures. Landslide area is distinguished from the non-landslide area by its various grades of materials, lower quartz content, higher clay content, finer composition of parental rocks and soils, higher average ratio of smectite versus quartz, lower compressive strength, lower shear strength, lower angle of internal friction, lower cohesion value, more poorly graded, higher soil moisture, higher liquid limit and plasticity index, and factor of safety is generally less than 1. The result of the sensitivity analysis indicated that in a failure condition, the cohesion and angle of internal friction angle values are equal to 2.51 kPa and 13.700, respectively. However, the results of laboratory tests showed that the cohesion and angle of internal friction values are equal to 8.59 kPa and 14.20. These results verified the landslide susceptibility map that the high landslide susceptibility area has the engineering characteristics of the landslide-prone area, while the low susceptible one has the characteristics of the non-landslide area 1. Introduction related to fundamental factors, such as soil and A landslide is the movement of a mass of rock, rocks properties, geological structures, and debris, or earth down a slope under the influence of topography, and the second group consists of gravity (Guzzetti, 2005). Landslide hazard is one of external factors such as anthropogenic and rainfall intensity. Kim and Song (2015) explained that the main natural disasters that need attention in Thailand. In 2006, two main landslide events caused landslides might occur when the external factors, as a trigger, overlap with the inherent factors. For a massive loss of both lives and properties in Uttaradit and Nan Provinces. When the Uttaradit example, heavy rainfall may trigger landslides on an unstable clayey soil slope. Heavy rainfall, in this Province landslide happened on 22-23 May 2006, the damaged cost was 308,615,331 Baht or roughly case, works as the external factor, while the unstable equal to 10 million US$ (30 Baht = 1 US$)(Usamah slope is the fundamental factor. The mentioned and Arambepola, 2013). The condition of this catastrophic translational landslides are shown in happened landslide in Uttaradit was mainly Figure 1. Laplae District is the most north-west influenced by the rainfall and the physical properties district in Uttaradit Province (as shown in the Figure of soil (Protong et al., 2018). The maximum rainfall 2), where it was also swept away by the 2006 intensity in May 2006 when the landslide occurred landslide. The lithology that composed of this study was 330.0 mm/day, and the limit for maximum area is Greywacke, Shale, Claystone, and mixture of rainfall intensity before the villagers need to shale and volcanic ashes (Khositanont et al., 2016 ). evacuate is 150.0 mm/day (Mairaing et al., 2016). Many previous studies (Fell et al., 2008, 1994, Similarly, like the previously mentioned case, the Hwang et al., 2009 and Kim and Song, 2015) mentioned that the influencing factors on landslides rainfall works as the external factor and the physical could be divided into two groups. The first group is properties of soil work as the inherent factor. Several studies showed that particular physical properties, such as grading, unit weight, liquid limit, International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

and plasticity indexes, are noted as factors affecting methodology in Figure 4. This study not only slope stability as well (Fonseca et al., 2017, Kim evaluates and compares the geotechnical properties and Song, 2015 and Mugagga et al., 2012). The and geology of landslide locations but also works as an engineering verification of a landslide status in physical properties of soil depend on the weathered the GIS-based landslide susceptibility mapping of grade and geology (Hurchinson, 1988, Lumb, 1975, the area (Herianto, 2020) between the past location Matsushi et al., 2006, Wakatsuki et al., 2005 and of landslides and landslide-free (or commonly called Yalcin, 2007). The characteristics of the weathered as non-landslide) location. Moreover, in order to soil are derived from the characteristics of its compare and to evaluate these parameters on parental rocks. landslides, the properties of the soils that cover slopes are measured, and the geological conditions Living in the landslide-prone area means that the are investigated. local people live in a high-risk area, especially when heavy downpours last for more than 24 hours like Three study areas are selected to consider and to what happened in 2006. As this natural disaster validate different statuses of their historical landslide in 2006 based on scars data and previous frequently occurs in Laplae District, the researcher landslide susceptibility map (as shown in geotechnical and geological conditions are the Figure 3). All of three chosen area (as shown in fundamental parameters that are scientifically and Figure 2 and 3) have different type of landslide socially important. Several researchers did their status (based on the scar data and the previous researcher’s landslide susceptibility map), which are research about the Laplae Landslide (Moazzam et a landslide area (very high susceptibility), a non- al., 2007, Phattaraporn et al., 2017, Protong et al., landslide area, and a non-landslide area close to the 2018, Tanang et al., 2010 and Usamah and landslide area (low susceptibility). A detailed Arambepola, 2013) but none have ever done about using geological and geotechnical characteristics in geotechnical and geological survey is performed as order to verify the landslide susceptibility mapping well as soil and rocks sampling in each study area. A of the area. Herianto (2020) mentioned that they used the series of laboratory soil tests are then conducted to frequency ratio method in order to make the measure the soil properties. Based on these results, landslide susceptibility map, which has a range accuracy of the map from 49 to 89%. However, this the relationships between landslide occurrence and accuracy verification was done by using area under soil properties are investigated and evaluated, and curve and receiver operating characteristics curve, a the landslide susceptibility map will be well verified statistical approach. Commonly, engineers prefer to after the field verification is done. Hutchison (2014) validate further by using field and laboratory agrees that Laplae is roughly located in the Sukhotai measurements to know whether the map that has Zone, where it is bounded with Nan-Uttaradit Suture been created is accurate enough based on physical on the East part which extends on the North-South field data. Herianto (2020) divided the landslide susceptibility map of Laplae area into five main direction, Chiang Mai Suture on the West part landslide susceptibility index (LSI), as shown in which also extends from North to South, Mae Ping Figure 3. The lower value of the LSI, meaning the Fault on the South part which azimuth is trending lower probably a landslide will occur, while the from North- West to South-East and Ailao Shan-Red higher the value of the LSI, meaning the higher probably a landslide will occur. Those five classes River Shear on the North part which also extends were subjectively renamed as very low from North-West to South-East. susceptibility, low susceptibility, medium susceptibility, high susceptibility, and very high Khositanont et al., (2016) divided Laplae area susceptibility. into three different formations, geologically. There are Lower Lap-Lae Formation (Pll1), Upper Lap-Lae However, the susceptibility maps formed was Formation (Pll2), and Alluvial Deposits (Qa). The verified using only existing landslide location (scar Lower Lap-Lae Formation consists of meta- data) by a statistical approach. In order to better understand and to verify the status of the landslide sandstone, greywacke, shale, and basaltic andesite susceptibility map in the field, direct field and sills. This Formation is called as Pll1. The Upper laboratory approaches will be done accordingly. The Lap-Lae Formation, on the other hand, consists of direct validation not only is important for increasing meta-shale or mudstone, greywacke sandstone, and the susceptibility map precision, but its data can also mudstone. This Formation is called as Pll2. be used for future landslide susceptibility modeling as input data, as shown in the flowchart of Geologically, those two Formations have been continuously weathered and eroded, which then become the source of the Quaternary deposits/soils. International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

The Alluvial Deposits consist of gravel, sand, silt, rainfall-triggered landslides. Huai Tong Sad (HTS) and clay. The age of this part is Quaternary, which is 0.01-1.6 Million years ago. is in the Center part of Mae Phun, Laplae District, Thailand, with elevation ranges between 100-400 m 1.1 Study Sites above sea level. Most of this area is used for Ban Dee (BD) is in the North-East part of Mae Phun, Laplae District, Thailand, with elevation harvesting, the variation of productive trees ranges between 200-500 m above sea level. Most of dominantly are durian trees and langsat trees. There this area is still on their natural state, and even there are also some parts of the area that are inhabited by are some cut roads established for local local villagers’ houses. This area was reported to be transportation, the vegetation dominated are natural the most impacted area by the 2006 rainfall- trees with large and deep roots. This area was triggered landslides. Not only be the most impacted reported to be the less impacted area by the 2006 area of landslide, but this area is also one of the most populated areas in Mae Phun. Figure 1: The catastrophic Laplae landslide (Mairaing et al., 2016) Figure 2: Research Location accessed in Google Earth International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Figure 3: Landslide susceptibility index map of Mae Phun, Uttaradit, Thailand (Herianto, 2020) Landslide Susceptibility Modelling Accuracy Assessment Statistical Validation Field-check Validation Is the map precise enough? Previous research No Yes Part of research Verified Landslide Susceptibility Map Figure 4: Methodology for landslide susceptibility mapping International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Saphan Huai Tong Sad (SHTS) is located 200-meter distribution curve is used to calculate the coefficient North of the Huai Thong Sad site, with elevation of uniformity and the coefficient of concavity. Not ranges between 100-400 m above sea level. Most of only indicating their engineering properties, but the this area is also being used for harvesting. This area, gradation of soil is also related to compressibility, however, was reported to be less impacted area by friction angle, and shear strength (Islam et al., 2011 the 2006 rainfall-triggered landslides. The anomaly and Mostefa et al., 2013). The size of particles is of being next to the most impacted area of landslide directly proportional to the peak friction angle and despite not significantly impacted by a landslide is angle of internal friction. The peak friction angle something that makes this area appealing to have a increases when the size of particles increases. more in-depth investigation. Atterberg limits are essential properties to point 2. Methodology out soil expansion potential at different moisture and clay contents (Selby, 1993). Those properties 2.1 Field Surveys explain the slope’s susceptibility to various slope Field surveys entailed mapping the Laplae area activities .The Atterberg’s limits were determined using software such as Rockd (Wisconsin-Madison, (ASTM, 2017b) .The higher clay content in the soil 2019) and SuperSurv (Inc., 2019). A total of three makes the higher PI in soil (Mugagga et al., 2012). sites was undertaken using a geological compass and Furthermore, the plasticity of the soils was further a tape measure. At each point, pits were dug to a finalized using the unified soil classification system depth of 100 cm. Undisturbed soil samples were (USCS) in the plasticity chart, which also enabled obtained at different depths (generally at an interval further classification of the fine material. The multi- of 30 cm) using KU-miniature sampler, while stage direct shear test uses a single soil specimen disturbed soil samples were obtained at different and shears the sample in stages with increasing depths of 30 cm also by using shovels. These normal stresses to minimize the testing sample samples were kept in air-tight-zipped plastic bags. (Mairaing et al., 2005). In addition, the multi-stage test is not an ASTM standard method for obtaining Undisturbed soil samples were collected by using total or effective stress parameters but has been KU- miniature thin-wall sampler, developed by widely used in practice (Nam et al., 2011). Kasetsart University (Mairaing et al., 2005), as Basically, the Mohr-Coulomb envelope is defined as shown in Figure 5a. These soil samples were used the relationship of strength parameter (c’, ø′). The for the determination of cohesion and angle of three samples were tested in the conventional internal friction under normal loading by conducting automatic direct shear test instrument, as shown in a multi-stage direct shear test in a consolidated Figure 5b and 5c. This study applied a multi-stage drained test. direct shear test in order to find the shear strength of undisturbed samples. It was developed from a direct Not only soil samples that were taken, the rocks shear test. It is said that a multi-stage direct shear samples and geological structures were also test is well suited for unsaturated soil testing where obtained and identified by using geological compass every single stage is tested until it is close to the by finding out strike/dip of the rock strata, rocks failure point in each normal load (at least three or composition, and structures. Several field tests were four normal loads) (Thongkhao et al., 2012). conducted, such as pocket penetrometer (ASTM, 2010) for finding compressive strength and pocket After getting the strength parameters, cohesion, shear vane (ASTM, 2019) to find the shear strength. and angle of internal friction, those values then used to calculate the Fs of slopes. If the Fs is > 1, it 2.2 Laboratory Analyzes means that the slope is inherently stable. However, In order to characterize and find out the physical it might fail by external factors causes. On the other properties of the slope materials in terms of its hand, if the Fs is < 1, it means that the slope is implications for slope stability, a series of analyzes naturally unstable (Berry and Reid, 1987). The Fs of were carried out at the Geotechnical Laboratory, an infinite slope is calculated by using Equation 1: Civil Engineering Department, Faculty of Engineering, Kasetsart University Kamphaeng Saen ������������ = [1 − ������������ . ������������ ] tan ������ + ������������ ������ 2������ 2������ Campus. The analyzes focused on the particle-size ������������ ������ tan ������ ������ sin distribution of soils, Atterberg limits, and shear Equation 1 strength. Petrography analysis was conducted in order to define the mineral composition of soil where ������������ = depth of landslide, ������ = cohesion, ������ = parental-rocks. Particle-size distribution test was slip depth, ������ = angle of internal friction, ������ = slope conducted (ASTM, 2017c). The particle-size angle, ������������ = unit weight of water (1g/cm3), ������������ = unit International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

weight of soil, ������ = gravitational acceleration (9.81 sorted, ¼ - 2 mm, consist of 50% K-feldspar, 20% m/s2). A back analysis is also done in order to Ca/Na-Feldspar, 10% micas, 5% quartz and carbonates compare the value of C and ������ from that with those cement of 15%, massive, good porosity. Total obtained in the laboratory. The effect of C and ������ on the factor of safety was also given by doing the thickness is more than 3 meters. The thickness of sensitivity analysis. the carbonaceous sandstone outcrop in Figure 8 is approximately 5 meters. The color of the soils of 3. Results and Discussion sandstone is brownish-white. The grain size particle 3.1 Field Surveys and Mapping is between <1/16 mm to 2 mm, well-sorted, sub- The results from the field observations showed that, rounded, grain-supported. Lithic fragments consist at the HTS landslide site, a translational debris flow of sandstone: whitish color, well-sorted, ¼ - 2 mm, had occurred. Evidence from past scars at the HTS consist of 60% Ca/Na-feldspar, 20% K-Feldspar, 10% site indicated an intermixture of materials between micas, 5% quartz and carbonates cement of 5%, deeper depth materials and shallower materials. The SHTS site, however, was the non-landslide location. massive, good porosity. Total thickness is more than A further remote sensing analysis was also 5 meters. accomplished in order to verify landslide scars location by using a 15 m DEM, satellite image, and Strikes and dips, as well as lineaments in the previous researchers’ mapping results (Herianto, location, are shown in Figure 9. These data are taken 2020 and Khositanont et al., 2016). The landslides were found to have occurred mostly on concave directly from the field. It is then combined with slopes ranging between 30-40 slope angles. The lithology, soils, and other topographic lineaments to average depth of landslide is 2.3 meters, and the make a geological map as well as to make a average depth of slip is 3 meters. The percentage of geological structure analysis. From the data, it landslides that occurred in the meta-shale unit is 35.44%, while 64.56% happened in the meta- shows that many lineaments of structures are greywacke unit. Field measurements revealed that controlling the research area. Most of the trend of the compressive strength and the shear strength of the lineaments are elongated in the North-west the BD site, HTS site, and SHTS site, respectively, South-east direction. From this data as well, faults are 539 kPa and 54 kPa; 294 kPa and 44 kPa; 490 and folds can be inferred based on the trend of kPa and 49 kPa. It agrees with other researchers’ structures, the trend of rivers direction, as well as results that the landslide area is having lower contours elevation difference. The data of compressive strength and shear strength than the geological structures data, as shown in Figure 9, is non-landslide area (Carruba annd Moraci, 1993, then plotted and analyzed in the stereo net and the Iannacchione and Vallejo, 2000 and Szafarczyk, rose-net, as shown in Figure 10 and Figure 11. It 2019). The location has elevation ranges from 100 shows from the rose-net analysis that the mean to 500 meters above sea level. The geomorphology direction of the force comes from the North-West of the sites is moderate-highly dissected direction, which then makes folding that have denudational hills. The main lithology of the North-East South-West folding plane direction. The mean direction of the fractures stress is N 342.60 E. outcrops in location mostly consists of chert, siltstone, and sandstone. The geological map and the geological section of the research location is shown in Figure 12 and The thickness of the outcrop in Figure 6 is Figure 13, respectively. The units of geology in the approximately 3 meters. The color of the soils from location are divided into three different units. The chert is dark brownish-red. The grain size particle is <1/16 mm, amorphous, layering structure, strike/dip oldest unit, which is part of Pll1, is consisted of N160E/550, mono-mineral silicas 100%, very low siltstone, in some places interbedded with sandstone and chert. porosity. Total thickness is more than 5 meters. The thickness of the carbonaceous siltstone outcrop in This unit is in green color on the map, which is Figure 7 is approximately 2 meters. The color of the called as Siltstone Interbedded with Sandstone and soils of siltstone is brownish-red. The grain size Chert Unit. The younger unit, which is part of Pll2, particle is between 1/256 mm to >64 mm, poorly is consisted of sandstone interbedded with chert. sorted, sub-angular, matrix-supported. Lithic fragments consist of sandstone: whitish color, well- This unit is in yellow color on the map, which is called as Sandstone Interbedded with Chert Unit. The last unit, which is part of Qa, is consisted of alluvial deposits. This unit is in orange color on the map, which is called as Alluvial Deposits Unit. International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

(A) (B) (C) Figure 5: (A) KU-miniature instrument (Mairaing, Thaijeamaree, & Kulsuwan, 2005), (B) the complete instruments of the direct shear test, and (C) the direct shear box equipment Figure 6: An Outcrop in BD area Figure 7: An Outcrop in HTS Area International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Figure 8: An Outcrop in SHTS Area Figure 9: Stike/Dip and lineaments map of laplae, Uttaradit, Thailand International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Figure 10: Stereo net analysis of bedding plane and geological structures Figure 11: Rose-net analysis of bedding plane and geological structures Figure 12: Geological map of laplae, Uttaradit, Thailand International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Figure 13: Geological section of laplae, Uttaradit, Thailand Figure 14: Shale thin section in the BD site under the normal-light microscope Figure 15: Siltstone thin section in the HTS site under the normal-light microscope International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

3.2 Petrography (Casini et al., 2011, Iannacchione and Vallejo, 2000, Petrography analysis has been done by using a Islam et al., 2011, Mostefa et al., 2013 and Stark et normal and cross-polarized microscope to show the al., 2014). Specifically, soils from BD site are mineral characteristics of soils’ parental rocks in greyish-brown gravelly sand, classified as well- each site with point counting or modal composition graded sand with little fines material (SW), while method. The analysis is only focusing on the mineral quartz, smectite, and illite as these three HTS soils are reddish-brown silty sand, classified as minerals are the most common minerals found in the poorly graded sand with little fines material (SP) parental rocks’ thin sections, while other minerals groups. The soils from SHTS are yellowish-brown such as Plagioclase, Lithic Fragment, Muscovite and Biotite are quantitatively counted, but it has a less gravelly sand, classified as well-graded sand with significant presence. little fines material (SW) groups. The average ratio of smectite vs. average of the 3.3 Soil Moisture and Atterberg Limits other two minerals: quartz and illite, in each site, Soil moisture and Atterberg limits, as shown in were also calculated. It shows that the average ratio Table 1, were measured and calculated for the of smectite with an average of quartz in HTS is behavior of soils in response to water content, and 37.10, while on the other sites, it falls below 1.3. the following implications for landslide occurrence. Moreover, the average ratio of smectite with an average of illite in HTS is 5.62, while on the other Even though the soils are coarse-grained material, sites, it only reaches 0.32, as shown in Table 1. It the Atterberg limits are calculated in order to means that smectite is the dominant minerals found analyze the behavior of the matrix-grained parts of to be the key mineral in the soil behavior in the the soil. The HTS soil has the highest soil moisture landslide area, while illite is the dominant minerals value. The liquid limit of HTS is above the in the non-landslide area. These results satisfy with threshold of 50% (ASTM, 2017a). It means that the other researchers’ result that smectite is found to be HTS clay has high swelling potential of the clay the key minerals in the landslide area (Ohlmacher, materials. The PI of HTS site also shows the highest 2001, Putra et al., 2019, Tamura and Hasegawa, 2015 and Zhao et al., 2007). expansion potential of the clay materials among other sites as well, it is considered as high plasticity The results show that the highest mineral clays. This means that the nature of the soils might composition percentage in each site is different. be one of the inherent factors. The plasticity chart, Smectite is found to be the most abundant mineral in HTS site, as shown in Figure 15, while in BD and which including PI that highlights the range of water SHTS, as shown in Figure 14 and Figure 16, Illite is contents where the soil exhibits plastic properties, is the most dominant mineral, as shown in Table 1.The shown in Figure 18. The plasticity chart shows that result was also verified by the Department of Mineral Resources of Kingdom of Thailand which even though all the soils of the non-landslide area is its laboratory analysis results show that the rocks critically near the 50% boundary of low and high found in BD, HTS, and SHTS sites are classified, plasticity, the only soil that contains high plasticity respectively, as shale, siltstone, and sandstone. clay is located at HTS site. This confirms that the 3.2 Soil Particle Distribution high plasticity of soils represents the characteristics Particle distribution curves are shown in Figure 17. of the landslide area as in accordance with others’ results (Bizimana and Sönmez, 2015, Kitutu et al., Soils from the respective sites were generally 2009 and Soralump, 2008). The nature of highly coarse-grained, with less than 50% of the material plastic soil of HTS site is displayed by their plotting passing the 0.075 mm sieve. D60, D30 and D10 of the above the boundary A-line and their PI is exceeded BD, HTS, and SHTS sites, respectively, are 5.50 the B-line threshold of 50%. mm, 2.00 mm, 0.40 mm; 2.00 mm, 0.80 mm, 0.40 mm; and 5.00 mm, 1.50 mm, 0.50 mm. Soils from 3.4 Shear Strength, Fs and Sensitivity Analysis The angle of internal friction and cohesion were the HTS site are finer and more poorly graded than computed by plotting shear strength versus normal those from the other two sites. The finer soil means stress curves, which were then used to calculate the Factor of Safety for each site. BD, HTS, and SHTS that the angle of internal friction and the peak sites’ curves are shown, respectively, in Figure friction angle is lower, which leads to a weaker 19(a,b,c). The number of points of scars vs. slope in strength of soil which characterized the landslide the research area is shown in Figure 20. The area, as confirmed with other researchers’ results landslides dominantly occurred in 300-400 slope, with a mean of 33.80. International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Figure 16: Sandstone thin section in the SHTS site under the normal-light microscope Table 1: Laboratory analysis results International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Figure 17: Particle distribution curve in the area Figure 18: Casagrande plasticity chart of the site location (Modified from Holtz and Kovacs (1981)) As expected by the Fs calculation of infinite slope slope angle of the area. When the slope angle of with various degree of slope, the slopes at BD and SHTS reaches 38, the Fs value falls below 1. It SHTS sites where the Fs > 1, tends to be more shows that even at the angle of 32 and 34, the slope stable, while the slopes at the HTS site where the at SHTS is still stable, but it becomes unstable when value of Fs is lower than the critical factor of 1, are it reaches an angle of 38 or more. The HTS landslide unstable. However, this only applies at the slope occurred on naturally unstable, inherently expansive angle of 32. The factor of safety also depends on the soils, harvested slopes. International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Figure 19: Shear stress versus normal stress curves for BD, HTS, and SHTS sites Table 2: Factor of safety and strength parameters at the three sites Dry density Slope Angle (0) Cohesion C Angle of Factor of (g/cm3) (kPa) internal safety 15.51 friction, ������ (0) Ban Dee 1.68 1.58 32 8.59 27.22 1.58 34 1.50 36 9.6 14.20 1.42 38 1.35 40 25.26 0.93 Huai Tong Sad 0.87 1.51 32 0.83 34 0.79 36 0.75 38 40 1.15 1.08 Saphan Huai Tong Sad 1.02 1.63 32 0.97 34 0.92 36 38 40 Number of points of scars Slope (degree) Figure 20: Frequency of landslide scars with various slope angle International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Figure 21: Sensitivity for slope Figure 22: Sensitivity graph of cohesion Figure 23: Sensitivity for friction angle It needs to be noted that even slopes where Fs > 1 is parameter varies, and other input parameters are kept constant in their mean values. The analysis conditionally stable, the external and internal factors assesses of which input parameter may be more critical to the assessment of the slope stability, and may change its stability when it applies their vice versa, which parameter has less effect on the influence on landslides such as rainfall intensity instability. The cohesion and friction angle of the failure surface was back analyzed by performing a rises, and the shrink-swell capacity of the clays sensitivity analysis. The effects of slope, cohesion, and angle of internal friction to the factor of safety triggered. These results also agree with other are presented in the form of sensitivity graphs in researchers’ results, which show that translational Figure 21, Figure 22, and Figure 23 in which the vertical axis represents the factor of safety and the failure mechanisms generally occurred at slopes horizontal axis represents the slope or cohesion or with angles between 80 to 320 (Schilter, 2019), angle of internal friction. between 300 to 450 (Huang et al., 2016), more than 250 (Zêzere, 2002). Sensitivity analysis assists researchers to evaluate the impact of an individual unknown variable, with the assumption that all other slope parameters are known. In this analysis, one International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Angle of Internal Friction (0) Cohesion vs Angle of Internal Friction FS 0.98-1.02 60 Friction= -4.7836 Cohesion + 59.041 50 R² = 0.9844 40 30 20 10 0 0 5 10 15 Cohesion (kPa) Figure 24: The relationship between cohesion and friction angle of the rock mass for an approximate factor of safety of 1, the dotted line is the line of factor of safety = 1 The value of cohesion and angle of internal friction materials, higher quartz content, lower clay content, by using the back analysis is compared with the more grainy composition of parental rocks and soils, value of cohesion and angle of internal friction lower average ratio of smectite vs. quartz, higher obtained from the laboratory.The slope, cohesion, compressive strength, higher shear strength, higher and the angle of internal friction obtained in the angle of internal friction, higher cohesion value, laboratory of the HTS site, where a landslide did more well-graded, lower soil moisture, lower liquid occur in 2006, is 300, 8.59 kPa and 14.20, limit and plasticity index, and factor of safety is respectively. On the other hand, from the back- generally more than 1. The factor of safety of SHTS analysis result, it is indicated that at the edge of area, however, falls below 1 once the slope angle failure, i.e., a factor of safety of 1, the slope, reaches 380. The result of the sensitivity analysis cohesion, and friction angle values were 29.630, indicated that in a failure condition, the cohesion 2.51 kPa and 13.700., respectively. The relationship and angle of internal friction angle values are equal of cohesion and angle of internal friction from the to 2.51 kPa and 13.700, respectively. However, the back analysis is shown in Figure 24. results of laboratory tests showed that in a failure condition, the cohesion and angle of internal friction 4. Conclusions values are equal to 8.59 kPa and 14.20. In a nutshell, The GIS-based landslide susceptibility map is now the hypothesis that geological and geotechnical field, and laboratory-verified that the high landslide properties as the inherent factors are related to the susceptibility area has the engineering landslide occurrence at the site where the slope characteristics of the landslide-prone area, while the failure occurred, is now generally accepted. low susceptible one has the characteristics of the non-landslide area. High susceptibility landslide Acknowledgments area is characterized by the intermixture of materials This research was funded by Kasetsart University of deep and shallow materials, lower quartz content, Postgraduate Scholarship Program for ASEAN higher clay content (especially smectite group clay), Students 2017, Kasetsart University Scholarships finer composition of parental rocks and soils, higher for ASEAN for Commemoration of the 60th average ratio of smectite vs. quartz, lower Birthday Anniversary of Professor Dr. Her Royal compressive strength, lower shear strength, lower Highness Princess Chulabhorn Mahidol, grant angle of internal friction, lower cohesion value, number 1/2560. This research was also supported by more poorly graded, higher soil moisture, higher the Department of Civil Engineering, Faculty of liquid limit and plasticity index, with the slope angle Engineering at Kamphaeng Saen, Kasetsart between 300 to 400 and factor of safety is generally University, Kamphaeng Saen Campus, Thailand. less than 1. Vice versa, low susceptibility landslide This research was also supported by a grant from area is characterized by more homogenous the Faculty of Engineering at Kamphaeng Saen, International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Kasetsart University, Kamphaeng Saen Campus, Fonseca, L., Luiz Lani, J., Fernandes Filho, E., Thailand. The author would like to thank the reviewers. Santos, G., Ferreira, W. and Maria Rocha Trancoso Santos, A. (2017). Variability in Soil Physical Properties in Landslide-Prone. Acta References Sci., Agron., Vol. 39, No.1. Kitutu, M. G., Muwanga, A., Poesen, J. and ASTM, 2010, New Test Method for Pocket Deckers, J. A., 2009, Influence of Soil Properties Penetrometer Test. In ASTM WK27337. West on Landslide Occurrences in Bududa District, Conshohocken, PA: ASTM International. Eastern Uganda. African Journal of Agricultural ASTM, 2017a, Standard Practice for Classification Research, Vol. 4, 611-620. of Soils for Engineering Purposes (Unified Soil Guzzetti, F., 2005, Landslide Hazard and Risk Classification System). In ASTM D2487-17. Assessment: Concepts, Methods and Tools for West Conshohocken, PA: ASTM International. the Detection and Mapping of Landslides, for ASTM, 2017b, Standard Test Methods for Liquid Landslide Susceptibility Zonation and Hazard Limit, Plastic Limit, and Plasticity Index of Assessment, and for Landslide Risk Evaluation. Soils. In ASTM D4318-17e1. West (ERLANGUNG DES DOKTORGRADS (DR. Conshohocken, PA: ASTM International. RER. NAT.)). RHEINISCHEN FRIEDRICH- ASTM, 2017c, Standard Test Methods for Particle- WILHELMS-UNIVESTITÄT BONN, Size Distribution (Gradation) of Soils Using Herianto, P., 2020, Landslide Susceptibility Sieve Analysis. In ASTM D6913 / D6913M-17. Mapping by Integrating Multidisciplinary Data. West Conshohocken, PA: ASTM International. (Master of Engineering). Kasetsart University ASTM, 2019, Standard Test Method for Kamphaeng Saen Campus, Thailand. Approximating the Shear Strength of Cohesive Holtz, R. D. and Kovacs, W. D., 1981, An Soils by the Handheld Vane Shear Device. In Introduction to Geotechnical Engineering: ASTM D8121 / D8121M-19. West Prentice-Hall. Presses Internationales Conshohocken, PA: ASTM International. Polytechnique. Berry, P. L. and Reid, D., 1987, An Introduction to Huang, W., Leong, E. and Rahardjo, H., 2016, Soil Mechanics: McGraw-Hill. Governing Failure Mode of Unsaturated Soil Bizimana, H. and Sönmez, O., 2015, Landslide Slopes Under Rainwater Infiltration. E3S Web of Occurrences in the Hilly Areas of Rwanda, their Conferences, 9, 15008. Causes and Protection Measures. Disaster doi:10.1051/e3sconf/20160915008. Science and Engineering, Vol. 1(1), 1-7. Hurchinson, J., 1988, General report: Carruba, P. and Moraci, N., 1993, Residual Strength Morphological and Geotechnical Parameters of Parameters from a Slope Instability. Proceedings Landslides in Relation to Geology and Third International Conference on Case Hydrology. Paper presented at the V Histories in Geotechnical Engineering, St. International Symposium on Landslides, Louis. 1481-186. Lausanne. 3-35. Casini, F., Brauchli, S., Herzog, R. and Springman, Hutchison, C. S., 2014, Tectonic Evolution of S. M., 2011, Grain Size Distribution and Particle Southeast Asia. Bulletin of the Geological Shape Effects on Shear Strength of Sand- Gravel Society of Malaysia, Vol. 60, 1-18. Mixtures. Paper presented at the 15th European Hwang, S., Guevarra, I. F. and Yu, B., 2009, Slope Conference on Soil Mechanics and Geotechnical Failure Prediction Using a Decision Tree: A Engineering (XV Case of Engineered Slopes in South Korea. ECSMGE).http://hdl.handle.net/20.500.11850/39 Engineering Geology, Vol. 104(1), 126-134. 588. doi:https://doi.org/10.1016/j.enggeo.2008.09.004 Fell, R., 1994, Stabilization of Soil and Rock Iannacchione, A. T. and Vallejo, L. E., 2000, Shear Slopes. Paper presented at the Proc. East Asia Strength Evaluation of Clay-Rock Mixtures. In Symp. and Field Workshop on Landslides and Slope Stability 2000. 209-223. Debris Flows, Seoul. Inc., S. T., 2019, SuperSurv (Version 10.0.0001) Fell, R., Corominas, J., Bonnard, C., Cascini, L., [iPhone software]: Supergeo Technologies Inc. Leroi, E. and Savage, W. Z., 2008, Guidelines Islam, M., Siddika, A., Hossain, M., Rahman, A. for Landslide Susceptibility, Hazard and Risk and Asad, M. A., 2011, Effect of Particle Size on Zoning for Land-Use Planning. Engineering the Shear Strength Behaviour of Sands. Geology, Vol. 102(3), 99-111. Australian Geomechanics Journal, Vol. 46, 85- doi:https://doi.org/10.1016/j.enggeo.2008.03.014 95. . International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Khositanont, S., Wattanapradda, S., Kunttasub, P., Nam, S., Gutierrez, M., Diplas, P. and Petrie, J., Suthisook, T., Dampitak, S., Nilda, K., 2011, Determination of the Shear Strength of Kaewmanee, P., and Pacca, I. (2016). The Study Unsaturated Soils Using the Multistage Direct On Geological Environments and Rainfall in Shear Test. Engineering Geology. Vol. 122(3-4), Relation to the Effect of Landslide. Research 272-280. and Development for Landslide Protection on Highland Slope Project Under the King Ohlmacher, C.G. 2001. The Relationship Between Initiation (Chaipattana Foundation). Geology and Landslide Hazards of Atchison, Kansas, and Vicinity. Current Research in Earth Kim, K.-S. and Song, Y.-S., 2015, Geometrical and Sciences, 244. Geotechnical Characteristics of Landslides in Korea Under Various Geological Conditions. Phattaraporn, S., Anusorn, R. and Thitawadee, S., Journal of Mountain Science, Vol. 12(5), 1267- 2017, Analyzing the Effects of Land Use changes 1280. for Landslide Susceptibility Assessment: A Case Study of LabLae district, Uttaradit province, Lumb, P., 1975, Slope Failures in Hong Kong. Thailand. Paper presented at the 38th Asian Quarterly Journal of Engineering Geology and Conference on Remote Sensing (ACRS 2017), Hydrogeology, Vol. 8(1), Delhi. doi:10.1144/GSL.QJEG.1975.008.01.02 Protong, S., Carling, P. A. and Leyland, J., 2018, Mairaing, W., Thaijeamaree, N. and Kulsuwan, B., Climate Change and Landslide Risk Assessment 2005, Landslide study in Phetchabun and in Uttaradit Province, Thailand. Engineering Chanthaburi. (Kasetsart University). Journal, Vol. 22(1), 243-255. doi:10.4186/ej.2018.22.1.243. Mairaing, W., Thaiyuenwong, S., Kulsuwan, B., Chantasorn, M., Pitayaniwit, J., Akejit, M., But- Schilter, J., 2019, Identifying Key Factors Affecting am, P., and Phumchaem, S. (2016). Geotechnical Translational Landslides in Part of the Yakima Engineering Modelling for Prediction of Fold and Thrust Belt, Washington State (Master Landslide Prone Area. Research and of Science (MS)). Central Washington Development for Landslide Protection on University, Highland Slope Project Under the King Initiation (Chaipattana Foundation). Selby, M. J., 1993, Hillslope Materials and Processes / M.J. Selby; with a Contribution by Matsushi, Y., Hattanji, T. and Matsukura, Y., 2006, A.P.W. Hodder. Oxford, England: Oxford Mechanisms of Shallow Landslides on Soil- University Press. Mantled Hillslopes with Permeable and Impermeable Bedrocks in the Boso Peninsula, Soralump, S., 2008, Landslide Hazard Investigation Japan. Geomorphology, Vol. 76(1), 92-108. Aand Evaluation of Doi Tung Palace: Example doi:https://doi.org/10.1016/j.geomorph.2005.10. of soil Creep Landslide. Paper Presented at the 003. EIT-JSCE Joint International Symposium \"Monitoring & Modelling in Geo-Engineering Moazzam, M. F. U., Vansarochana, A., Bangkok, Thailand. Boonyanuphap, J. and Choosumrong, S., 2017, Landslide Assessment Using Gis-Based Stark, N., Hay, A. E., Cheel, R. and Lake, C. B., Frequency Ratio Method: A Case Study of Mae- 2014, The Impact of Particle Shape on the Angle Phun Sub-District, Laplae District, Uttaradit of Internal Friction and the Implications for Province, Thailand. Paper presented at the 38th Sediment Dynamics at a Steep, Mixed Sand- Asian Conference on Remote Sensing (ACRS Gravel Beach. Earth Surf. Dynam., Vol. 2(2), 2017), Delhi. 469-480. doi:10.5194/esurf-2-469-2014. Mostefa Kara, E., Meghachou, M. and Aboubekr, Szafarczyk, A., 2019, Stages of Geological N., 2013, Contribution of Particles Size Ranges Documentation on the Example of Landslides to Sand Friction. ETASR - Engineering, Located on the Slopes of the Dam Reservoir Technology & Applied Science Research, Vol. “Swinna Poreba” (Poland). IOP Conference 3(4), 5. Series: Earth and Environmental Science, 221, 012037. doi:10.1088/1755-1315/221/1/012037. Mugagga, F., Kakembo, V. and Buyinza, M., 2012, A Characterisation of the Physical Properties of Tamura, E. and S. Hasegawa (2015). Verification of Soil and the Implications for Landslide swelling and landslides of smectite-bearing Occurrence on the Slopes of Mount Elgon, ground due to hydrothermal alteration in non- Eastern Uganda. Natural Hazards, Vol. 60(3), volcanic region. 10th Asian Regional 1113-1131. doi:10.1007/s11069-011-9896-3. Conference of IAEG. Kyoto, Japan. Tanang, S., Sarapirome, S. and Plaiklang, S., 2010, Landslide Susceptibility Map of Namli Watershed, Uttaradit, Thailand. Paper presented International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

at the 31st Asian Conference on Remote Sensing Wisconsin-Madison, U., 2019, Rockd (Version 2010. 2.6.0) [iPhone software]: UW Board of Regents, Thongkhao, T., Choowong, M., Charusiri, P. and 2019. Thitimakorn, T., 2012, Engineering Properties of Residual Soil from Landslide Hazard Area in Yalcin, A., 2007, The Effects of Clay on Landslides: Ban Nong Pla Village, Chiang Klang District, A Case Study. Applied Clay Science, Vol. 38(1), Nan Province, Northern Thailand. Bulletin of 77-85. Earth Sciences of Thailand. doi:https://doi.org/10.1016/j.clay.2007.01.007. Usamah, M. and Arambepola, N., 2013, Lessons Learned from the 2006 Flashfloods and Zhao, Y., et al. 2007. Evolution of Shear Strength of Landslide in Uttaradit and Sukhothai Provinces: Clayey Soils in a Landslide Due to Acid Rain: A Implication for Effective Landslide Disaster Case Study in The Area of Three Gorges, China. Risk Management in Thailand. Landslide First North American Landslide Conference, Science and Practice, Vol. 6. doi:10.1007/978-3- Landslides and Society: Integrated Science, 642-31319-6_88. Engineering, Management, and Mitigation. Vail, Wakatsuki, T., Tanaka, Y. and Matsukura, Y., 2005, Colorado. Soil Slips On Weathering-Limited Slopes Underlain by Coarse-Grained Granite or Fine- Zêzere, J. L., 2002, Landslide Susceptibility Grained Gneiss Near Seoul, Republic of Korea. Assessment Considering Landslide Typology. A CATENA, Vol. 60(2), 181-203. Case Study in the Area North of Lisbon doi:https://doi.org/10.1016/j.catena.2004.11.003. (Portugal). Nat. Hazards Earth Syst. Sci., Vol. 2(1/2), 73-82. doi:10.5194/nhess-2-73-2002. International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Step-wise Overlay Technique for the Mapping of Unconfined Groundwater Potential Zone in Tectonically Controlled Landforms Ihsan, H, M.,1 Hadi, M, P.2 and Sartohadi, J.3 1Graduate School of Universitas Gadjah Mada, Indonesia 2Faculty of Geography, Universitas Gadjah Mada, Indonesia 3Faculty of Agriculture, Universitas Gadjah Mada, Indonesia Abstract The existence of unconfined groundwater in tectonically-controlled volcanic landforms is continuous and discontinuous as in alluvial landforms. DEM-NAS is a product available all over Indonesia and has the potential for unconfined groundwater mapping. The study was conducted by applying GIS techniques to obtain information on drainage density, TWI (Topographic Wetness Index), slope angle, and lineament density with openly available algorithms. All parameters are put on the map using fuzzy analysis techniques before being combined using a step-wise overlay technique. Classification of unconfined groundwater potential was done based on the value of the merging parameter map. The results of the field verification show that 44 springs are in the potential zone with an area percentage of 68.01%, and 10 springs are in the less potential zone with an area percentage of 31.99%. 1. Introduction The existence of unconfined groundwater in Elevation Model) is a product that can be processed tectonically-controlled volcanic landforms is into parameters for the existence of potential continuous and discontiouos as in alluvial groundwater with the help of spatial technology landforms. Medium to low groundwater potential is (Elbeih, 2015). DEM can obtain terrain aspects usually found in weathered and cracked hard rocky such as drainage density, slope angle, lineament areas with high elevation and very steep density (Al-Shabeeb et al., 2018) and TWI topography, while flood plain zones are usually (Topographic Wetnes Index) (Setiawan et al., associated with good groundwater potential due to 2019). high infiltration levels of alluvium sediment (Thapa The study of groundwater potential using the et al., 2017). Volcanic landforms have high step-wise overlay approach (tiered overlapping elevation and steep topography. The tectonical analysis) has never been done in a tectonically- activity in volcanic landforms causes hard material controlled area. The mapping method of to slowly weathered. The existence of groundwater groundwater potential is usually performed by storage is usually in volcanoes associated with several methods, such as: the MCDA (Multi- sediment, especially pumice (Singh et al., 2019). Criteria Decision Analysis) method carried out in The existence of groundwater is affected by Guna tana area in the upper blue Nile Basin, geomorphological units such as plains, terraces, Ethiopia (Andualem and Demeke, 2019), the hybrid volcanoes, hills, valleys and material conditions computational intelligence model carried out in (Sandoval and Tiburan, 2019). Vadora District, Gujarat, India (Pham et al., 2019), The study of groundwater can be identified by GA (Genetic Algorithm) method performed in hydrological, geological, geophysical, and remote Wuqi Country, Shaanxi Province, China, and sensing approaches (Venkateswaran and comparisons of groundwater potential methods Ayyandurai, 2015). Exploration of groundwater between AHP (Analytic Hierarchy Process), through hydrogeological, geological and Catastrophe and entrophy techniques in geophysical techniques requires a great deal of cost Trichinopoly, India (Jenifer and Jha, 2017). The and time (Jha et al., 2010, Razandi et al., 2015 and whole method is carried out by utilizing remote Murmu et al., 2019). Remote sensing and GIS sensing technology and GIS. (Geographic Information System) are efficient tools The objective of this study was to apply the in mapping groundwater potential by reducing time step-wise overlay technique in the tectonically- and energy (Nithya et al., 2019). DEM (Digital controlled landforms in Garut Basin, West Java International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Province, Indonesia for mapping the unconfined contact, so this area is a tectonically-controlled groundwater potential. The landforms studied are in volcano quarter (Figure 2). the form of hilly zones surrounded by young volcanoes so that the spatial distribution patterns The study area is divided into 3 areas (Table 1), are distinctly different. Tectonical activity is namely structure 1, structure 2 and structure 3 characterized by the existence of faults, touches, (Figure 1). Structure 1 is the landform dominated synclines and anticlines in the study area. The study by highlands but there is a rain shadow area with an parameters used were drainage density, TWI, slope area of 129 km². Structure 2 is the landform angle, and lineament density. The DEM used is dominated by lowlands, this area is adjacent to the DEM-NAS, available for data studies throughout landforms of volcanic plains with an area of 57.66 Indonesia. This study is expected to be useful for km ². Structure 3 is the landform dominated by the effective and efficient study of groundwater semi-highlands but it is not exist rain shadow area potential in similar landforms. with an area of 42.78 km². 2. Study Area 3. Data and Method The area of tectonically-controlled landforms is in The study of groundwater potential using DEM the northwestern part of the Garut Basin with an data produced by BIG (Geospatial Geoinformation area of 229.54 km². Garut Basin is located in Garut Agency) is accessed on http://tides.big.go.id- Regency, West Java Province, Indonesia. The /DEMNAS/Jawa.php and available for all areas in Garut Basin is surrounded by a volcanic cone Indonesia. BIG stated that the DEM is called DEM- morphology. These conditions cause the Garut NAS which has a spatial resolution of 0.27 Basin to have a varied topography, slope angle, and arcsecond. DEM-NAS is used to generate elevation, from flat to mountainous. The basic groundwater parameters such as river flow density, materials of this area are andesite and basalt. Based TWI (Topographic Wetness Index), slope angle, on the Garut-Pamempeuk regional geological map, and lineament density.Data developed for zoning there are faults, synclines, anticlines and geologic analysis of groundwater potential are morphology, lithology and hydrology (Setiawan et al., 2019). Figure 1: Study site of tectonically-controlled landforms Table 1: Characteristics of the study area Location Elevation Morphology Relief Existing Structure 1 2238-746 Structural fault with a lava Dominated by Rain Shadow Sturcture 2 1648-596 dome Mountainous area (>40°) Area Structure 3 1062-613 Structural faults with Dominated by Hilly area Not a Shadow parasitic cones (15°-40°) Area Structural faults with Dominated by Wavy Not a Shadow Alluvial area (8°-15°) Area International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

DEM is able to generate hydrogeological data such Variable α (m²) and β (m) are performed based on as drainage density, slope angle (Das, 2019), TWI DEM extraction using SAGA GIS 2.3.2 software. (Topographic Wetness Index) (Misi et al., 2018), DEM is able to be used to develop TWI analysis and lineament density (Wilson and Gallant, 2000). (Sorensen and Sorensen, 2007). 3.1Drainage Density 3.3 Slope Angle Slope angle is one of the parameters that has a role Drainage density is the total length divided by the in recharge control of groundwater potential. Groundwater potential is usually at a low slope total area in each river flow (Horton, 1932). The angle value (Hammouri and Al-amoush, 2014). Steep slope angles produce small absorption, closeness of the channel and the nature of the because water flows in the direction of gravity according to its degree, so it does not have adequate surface material can be shown through the river time to infiltrate the surface and replenish the saturation zone (Yeh et al., 2016). The slope angle flow density in km/km² (Venkateswaran and is the main indicator of groundwater infiltration below the surface, usually on flat slopes and slow Ayyandurai, 2015). Areas with high river flow surface runoff, causing a lot of time for rain water to infiltrate and high infiltration power (Murasingh density values cause infiltration power to decrease et al., 2018). The slope angle can be generated through DEM data by extracting processed GIS and runoff to rise, so the groundwater potential is at software slope (Patra et al., 2018). The characteristics of the slope angle are described a low value of river flow density (Al-adamatand through the classification of Van Zuidam (1985) to clarify the description of the slope angle class in the Al-shabeeb, 2017). DEM is able to extract river study area. flow data with the direction algorithm method on the plugin arc hydro tool (Yousif et al., 2018). The results of river flow extraction will produce several river orders. The flow used is an order which is a perennial not periodic river. Drainage density can be calculated based on the following Horton formula (1932): ������ d = ∑ ������������/������(������������2) ������=1 Equation 1 3.4 Lineament Density The formula explains that d is the river density In the remote sensing description, lineament is often index (km/km²), Di is the length of the river including its tributaries (km²) and A is the area. reflected as ridge line and valleys with the River flow density can be analyzed using ArcGis software with the plugin arc hydro tools and traditional method of line extraction, which is based lineament density (Oikonomidis et al., 2015). on artificial visual or semi-automatic interpretations (Han et al., 2018). The lineament in this study is a straight-shaped feature associated with the valley, because the groundwater potential is in a flat area. 3.2 TWI (Topographic Wetness Index) The greater the value of lineament density or the The TWI study funtions to measure topographic control in the hydrological process (Setiawan et al., closer it is to large faults, the higher the 2019). TWI is able to describe the spatial pattern of water-saturated areas, which is the key to permeability (Al-Rozouq et al., 2019). GIS understanding the diversity of materials and hydrological processes (Grabs et al., 2009). The technology is able to identify the lines by utilizing higher the TWI value, the lower the slope angle value. and the more the groundwater potential DEM data (Rajasekhar et al., 2018). Lineament (Nejad et al., 2017). The physical properties of the soil and parent material have a correlation to the density can be calculated using the following TWI calculation results (Gillin et al., 2015). TWI calculationis carried out based on formulas made formula: by Beven and Kirby (1979): ������ = ∑������������==���1��� ������������ ������ Equation 3 TWI = In ( α ) Where LD is the lineament density∑������������==1������ ������������, is the total length of the lineament (L), and A is the area of tanβ landforms. High Ld value indicates the high porosity and groundwater potential (Yeh et al., Equation 2 2016). The lines can be analyzed using PCI Geomatica software. The lineament density is The formula explains that α is the area of the mapped using ArcGIS Software by utilizing the unslope contributing area and β is the slope angle. density function tools (Al-adamat and Al-shabeeb, 2017). International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Figure 2: Geology map of study area Figure 3: Drainage density Figure 4: Topographic wetness index Figure 5: Slope angle Figure 6: Lineament Density International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

3.5 Zoning of Groundwater Potential 1 for order 7-3 is 242.64 km, structure 2 for order 7- The method used for zoning groundwater potential 3 is 112.98 km and structure 3 for order 6-3 is 88.64 is stepwise overlay. This method is basically an km. The total area determines the length of the river overlay, but it is done in stages from one to four in tectonically-controlled landforms. The wider the parameters (Table 2). The overlay method in GIS area, the longer the river flows. induces complex matrix changes, produces binary changes in the image without changing the major Drainage density in tectonically-controlled minor landscape, and is based on classification landforms has a great value located in highlands, but (Tiede, 2014). Each tiered parameter overlay will springs are rarely found in structure 3 with a density show different results, but the more parameters, the value of 0 - 9.94 km². The structure 3 is has the better the accuracy assumed. lowest area and length of river, compared to other structures. The tectonical activity causes differences Table 2: Step-wise overlay of groundwater in the length and density of river flows, because the Parameters material in this area is brittle and it forms streams and river flows. The structure 1 has a river flow No Landforms Study Paramaters density value of 0 - 5.56 km², while the structure 2 is 0 - 8.99 km². Presentation of drainage density in Structure 1 spatial form illustrates the high and low of an area, so that the distribution of groundwater potential will 1 Structure 2 Drainage density (2) + TWI (1) appear based on parameters of drainage density (Figure 3) Structure 3 4.2 TWI (Topographic Wetness Index) 2 Structure 1 (Drainage density + TWI) (2) + The TWI study on tectonically-controlled landforms Structure 2 Slope (1) has a value with a small difference, even almost the Structure 3 same between structure 1, structure 2 and structure 3. In tectonically-controlled landforms, the TWI Structure 1 (Drainage density + TWI + value of structure 1 is 6.74-14.51, structure 2 is 6.68 - 13.68 and structure 3 is 6.72 - 14.29. This 3 Structure 2 Slope) (2) + Lineament density condition is caused by the formation process of the same land and material originating from volcanoes. Structure 3 (1) Based on the geological map of the Pamempeuk- Garut sheet from the Geological Agency, typical The processing results of each parameter of volcanic rocks such as andesite-basalt, breccias and unconfined groundwater (Drainage Density, TWI, tuffs are found. slope angle and lineament density) were raster data that will be processed through fuzzy membership to Surface material is produced from extrusive equalize raster values with a range of value 0-1 in volcanoes that surround the Garut Basin. The ArcGis software. Parameter weights were striking difference in TWI value depends on relief determined by trying to overlay each parameter in units such as upper volcanic slope, middle volcanic rotation until finding a good value (Table 2). slope, lower volcanic slope, caldera, old volcanic hills, old volcanic plains, old volcanic mountains, Zoning verification of groundwater potential valley floor, collovial slopes and subresent lava was done by surveying the location of a spring, then flows (Figure 4). the coordinates location was plotted. The sampling technique used was purposive sampling because the High TWI value is found in areas with flat existence of groundwater was unknown in number topography. A high TWI value is assumed to have a and location, so the numbers were determined based high infiltration rate, due to the existence of non- on field conditions. Zoning classification of steep terrain, allowing small runoff. The TWI study groundwater potential was divided into two, namely functions to measure topographic control in the potential area and less potential area. The step-wise hydrological process (Setiawan et al., 2019). TWI is overlay method was applied to each area of structure able to describe the diversity of materials and 1, structure 2 and structure 3. hydrological processes through the spatial pattern of water-saturated areas (Grabs et al., 2009). Based on 4.1 Drainage Density a 1: 50,000 semi-detailed soil map from BBLSDP Analysis of river flow using hydro arc tools shows (Center for Agricultural Research and Development that the shape of structure 1 and structure 3 has 7 on Land Resources), haplic andosol, district river orders while structure 2 has 6 river orders. In cambisol, haplic podzolic, rodic latosol and gleisol each form, the river order 1 and 2 are not used, district were found. The soil unit in the study area because it is assumed that the river is perennial. The level of infiltration in the area around the perennial river is greater than in the intermittent river. Groundwater potential is at a low drainage density level, due to low infiltration rates (Al-adamat and Al-shabeeb, 2017). The length of the river structure International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

has good drainage and smooth texture. lineaments and the area. Based on the analysis, 4.3 Slope Angle structure 1 has a lineament density of 0-6.75 km², In general, the slope angle in the tectonically- structure 2 has a lineament density of 0-6.62 km², controlled landforms is dominated by hilly to and structure 3 has a lineament density of 0 - 9.17 mountain relief classes (16°-35°). Classification of km². Areas with high lineament values are areas slope angles from hilly to mountainous (16°-35°) in with high groundwater potential, due to high structure 1 has an area of 60.89 km² and structure 3 permeability (Al-Rozouq et al., 2019). The has an area of 17.03 km². Based on the Van Zuidam distribution of lineaments density is depicted slope angle classification (1985) the area of through raster data which is an infographic of the structure 1 and structure 3 has 7 classifications, such existence of groundwater potential (Figure 6). as flat to almost flat (0°-2°), wavy (2°-4°), wavy to corrugated (4°-8°), corrugated to hilly (8°-16°), hilly 4.5 Unconfined Groundwater Potentials to mountainous (16°-35°), steep mountains (35°- DEM-NAS is able to analyze the study of 55°) and very steep mountains (>55°). The structure groundwater potential using parameters of drainage 2 has 6 slope angle classifications. density, TWI (Topographic Wetness Index), slope angle and lineament density. Based on field The classification of very steep mountain slope verification there are 54 springs in tectonically- angle (> 55°) is not found in the structure area 2. controlled landforms (Figure 7). Step wise overlay This condition is due to the fact that structure 2 is analysis shows 44 springs in potential landforms dominated by corrugated to hilly slope angle classes with an area of 156 km² (68.01%) and 10 springs in (8°-16°) with an area of 18.12 km². The structure 2 areas with less potential landforms with an area of is adjacent to the volcanic landforms. The condition 73.42 km² (31.99%) (Table 3). of the slope angle has groundwater potential because it has a longer residence time for rainwater The stepwise overlay method shows that and a high infiltration rate (Patra et al., 2018). The everytime each parameter is added, the number of highest slope angle is in the structure 1 which is springs in a potential landform increases or remain 59.25 ° (Figure 5). the same as the number of springs in potential landforms, while the total area increases or 4.4 Lineament Density decreases.The results of the step wise overlay The study of lineament applied PCI Geomatic analysis on stage 1 analysis (drainage density + software shows straight-shaped features on the TWI) has the same number of springs as stage 2 surface of the earth. This lineament is valley and analysis (drainage density + TWI + slope angle), ridge. In this study, the lineament used is valley which is 41 springs in potential landforms and 13 because the groundwater potential is in a flat area springs in less potential landforms (Table 2). This such as valleys. The extraction results were condition is because TWI and slope angles both manually sorted using ArcGIS based on the have topographic control functions of rainwater hillshade and color appearance of the DEM-NAS. A runoff and infiltration. low DEM-NAS value indicates low elevation or valley, while a high DEM-NAS value is the ridge. Stage 3 analysis is the result of step wise overlay Based on the analysis, 650 lineaments are found in for all parameters. Structure 1 has the highest structure 1 with a length of 301.28 km, 161 number of springs in the potential landform. This lineaments are found in structure 2 with a length of area is dominated by highlands and includes 67.76 km, and 170 lineaments are found in structure mountainous areas. In the western part of the 3 with a length of 75.85 km. The number of structure 1 there is a rain shadow area. Groundwater lineaments is affected by terrain conditions which sources are indicated to originate from orographic are undulating terrain and total area. rain that occurs in mountains and tectonically- controlled landforms. The lineament is rarely found in flat areas, due to less varied topographic conditions such as those Tectonism control is characterized by hot in the structure 2. The lineament can be an springs and waterfalls in the structure 1. Stage 1 indication of the existence of structure in an area. analysis (drainage density + TWI), stage 2 analysis The existence of structures must be further verified (drainage density + TWI + slope angle), and stage 3 based on field data and geological maps. Based on analysis (drainage density + TWI + slope angle + the geological map of the Pamempeuk-Garut sheet lineament density) show 20 springs in potential from the Geological Agency, the tectonically- landforms and 3 springs in less potential landforms controlled landform has a structure of faults, (Figure 6). Structure 2 and structure 3 experience anticline, syncline and geologic contact. The the addition of springs in the potential area after lineament density is determined by the number of stage 3 analysis. International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Table 3: Parameters of unconfined groundwater potential Spring Number of springs Landforms Study Parameters Potential Area Less Area (km²) 23 Structure 1 drainage density + TWI 20 (km²) Potential 33.46 23 Structure 2 16 18.03 8 Structure 3 Total 5 95.58 3 17.12 54 drainage density + TWI + 41 68.61 23 Structure 1 20 39.72 7 37.06 23 Structure 2 slope 16 8 Structure 3 Total 5 25.63 3 20.36 54 drainage density + TWI + 41 16.37 23 Structure 1 slope + Lineamen Density 20 160.93 13 73.79 23 Structure 2 Total 18 36.04 Structure 3 91.82 3 8 6 20.71 54 44 37.45 7 16.67 26.48 3 155.75 13 73.42 92.45 3 36.96 5 26.71 2 156.12 10 Figure 7: Distribution of unconfined groundwater potential zone International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Analysis on stage 1 and 2 did not experience a Ethiopia. Journal of Hydrology: Regional change in the number of springs in the potential and Studies, Vol. 24, 100610. less potential landforms (Table 2). Potential and less Beven, K. J. and Kirkby, M. J., 1979, A Physically potential areas at each stage of analysis differ, due Based, Variable Contributing Area Model of to the impact of the step-wise overlay results on Basin Hydrology. Hydrological Sciences raster data. Bulletin, Vol. 24, 43-69. Das, S., 2019, Comparison among Influencing 5. Conclusion Factor, Frequency Ratio, and Analytical DEM-NAS is able to analyze the study of Hierarchy Process Techniques for Groundwater groundwater potential using parameters of drainage Potential Zonation in Vaitarna Basin, density, TWI (Topographic Wetness Index), slope Mahasashtra, India. Groundwater for angle and lineament density in tectonically- Sustainable Development, Vol. 8, 617-629. controlled landforms. Step wise overlay analysis Elbeih, S. F., 2015, An Overview of Integrated shows that 44 springs are in the potential landforms Remote Sensing and GIS for Groundwater with an area percentage (68.01%) and 10 springs are Mapping in Egypt. Ain Shams Engineering in less potential landforms with an area percentage Journal, Vol. 6, 1-15. (31.99%). The greatest value of the river flow Gillin, C. P., Bailey, S. W., Mcguire, K. J. and density is located in the area of highland dominance Gannon, J. P., 2015, Mapping of but springs are rarely found, which is in structure 3. Hydropedologic Spatial Patterns in a Steep The condition of the relief unit and the process of Headwater Catchment. Soil Science Society of landform formation affect the TWI value America Journal, Vol. 79, 440-453. distribution. In general, the slope angle in the Grabs, T., Seibert, J., Bishop, K. and Laudon, H., tectonically-controlled landforms is dominated by 2009, Modeling Spatial Patterns of Saturated hilly to mountainous relief classes (16°-35°). The Areas: A Comparison of the Topographic lineament density is determined by the number of Wetness Index and a Dynamic Distributed lineaments and the area. The lineament of the Model. Journal of Hydrology, Vol. 373(1-2), 15- extraction results using spatial technology is a 23. valley or ridge in the tectonically-controlled Hammouri, N. and Al-amoush, H., 2014, landforms. Groundwater Recharge Zones Mapping Using GIS: A Case Study in Southern Part of Jordan References Valley, Jordan. Arab J Geosci, Vol. 7, 2815- 2829. Al-adamat, R. and Al-shabeeb, A. A., 2017, A Han, L., Liu, Z., Ning, Y. and Zhao, Z., 2018, Simplified Method for the Assessment of Extraction and Analysis of Geological Groundwater Vulnerability to Contamination. Lineaments Combining a DEM and Remote Journal Water Resource and Protection. Vol. 9, Sensing Images from the Northern Baoji Loess 305-321. Area. Advances in Space Research, Vol. 62, 2480-2493. Al-rozouq, R., Shanableh, A., Merabtene, T., Horton, R. E., 1932, Drainage Basin Characteristic. Siddique, M., Khalil, M. A., Idris, A. and Trans Am Geophys Union, Vol. 13, 350-361. Amulla, E., 2019, Potential Groundwater Zone Jenifer, M, A. and Jha, M. K., 2017, Comparison of Mapping Based on Geo-Hydrological Analytic Hierarchy Process, Catastrophe and Considerations and Multi-Criteria Spatial Entropy Techniques for Evaluating Groundwater Analysis: North UAE. Catena, Vol. 173, 511- Prospect of Hard-Rock Aquifer Systems. 524. Journal of hydrology, Vol. 548, 605-624. Jha, M. K., Chowdary, V. M. and Chowdhury, A., Al-shabeeb, A. A., Al-adamat, R., Al-fugara, A. and 2010, Groundwater Assessment in Salboni Al-amoush, H., 2018, Groundwater for Block, West Bengal (India) Using Remote Sustainable Development Delineating Sensing, Geographical Information System and Groundwater Potential Zones within the Azraq Multi-Criteria Decision Analysis Techniques. Basin of Central Jordan Using Multi-Criteria Hydrogeol. J. Vol. 18, 1713–1728. GIS analysis. Groundwater for Sustainable Misi, A., Gumindoga, W. and Hoko, Z., 2018, An Development, Vol. 7, 82-90. Assessment of Groundwater Potential and Vulnerability in the Upper Manyame Sub- Andualem, T. G. and Demeke, G. G., 2019, Catchment of Zimbabwe. Physics and Chemistry Groundwater Potential Assessment Using GIS and Remote Sensing: A Case Study of Guna Tana Landscape, Upper blueNile Basin, International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

of the Earth, Vol. 105, 72-83. Potential Mapping Using GIS. Earth Sci. India, Vol. 8, 867-883. Murmu, P., Kumar, M., Lal, D., Sonker, I. and Sandoval, J. A. and Tiburan, C. L., 2019, Singh, S. K., 2019, Deleniation of Groundwater Identification of Potential Artificial Potential Zones Using Geospatial Techniques Groundwater Recharge Sites in Mount Makiling and Analytichal Hierarchy Process in Dumka Forest Reserve, Philippines Using GIS and District, Jharkhand, India. Groundwater Analytical Hierarchy Process. Applied Sustainable Development, Vol. 9, 100239. Geography, Vol. 105, 73-85. Setiawan, O., Satohardi, J., Hadi, M. P. and Murasingh, S., Jha, R. and Adamala, S., 2018, Mardiatno, D., 2019, Delineating Spring Geospatial Technique for Delineation of Recharge Areas Infered from Morphological, Groundwater Potential Zones in Mine and Dense Lithological, and Hydrological Datasets on Forest Area Using Weighted Index Overlay Quaternary Volcanic Landscapes at the Southern Technique. Groundwater for Sustainable Flank of Rinjani Volcano, Lombok Island, Development, Vol. 7, 387-399. Indonesia. Acta Geophysica, Vol. 67, 177-190. Singh, S. K., Zeddies, M., Shankar, U. and Grif, G. Nejad, S. G., Falah, F., Daneshfar, M., Haghizadeh, A., 2019, Potential Groundwater Recharge A. and Rahmati, O., 2017, Deleniation of Zones within New Zealand. Geoscience Groundwater Potential Zones Using Remote Frontier, Vol. 10, 1065-1072. Sensing and GIS-based Data-Driven Models. Sorensen, R. and Seiber, J., 2007, Effects of DEM Geocarto International, Vol. 32, 167-187. Resolution on the Calculation of Topographical Indices: TWI and its Components. Journal of Nithya, C. N., Srinivas, Y., Magesh, N. S. and Hydrology, Vol. 347, 79-89. Kaliraj, S., 2019, Assessment of Groundwater Thapa, R., Gupta, S., Guin, S. and Kaur, S., 2017, Potential Zones in Chittar basin, Southern India Assessment of Groundwater Potential Zones Using GIS Based AHP Technique. Remote Using Multi-Influencing Factor (MIF) and GIS: Sensing Applications: Society. Environment, A Case Study from Birbhum District, West Vol. 15, 100248. Bengal. Appl. Water Sci., Vol. 7, 4117-4131. Tiede, D., 2014, A New Geospatial Overlay Method Oikonomidis, D., Dimogianni, S., Kazakis, N. and for the Analysis and Visualization of Spatial Voudouris, K., 2015, A GIS/Remote Sensing- Change Patterns Using Object-Oriented Data based Methodology for Groundwater Modeling Concepts. Cartography Geographic Potentiality Assessment in Tirnavos Area, Information Science, Vol. 41, 227-234. Greece. Journal of Hydrology, Vol. 525, 197- Van Zuidam, R. A., 1985, Aerial Photo- 208. Interpretation Terrain Analysis and Geomorphology Mapping. (The Hague: Patra, S., Mishra, P. and Chandra, S., 2018, SmithPublisher). Delineation of Groundwater Potential Zone for Venkateswaran, S. and Ayyandurai, R., 2015, Sustainable Development: A Case Study from Groundwater Potential Zoning in Upper Gadilam Ganga Alluvial Plain Covering Hooghly District River Basin Tamil Nadu. Aquatic Procedia, Vol. of India Using Remote Sensing, Geographic 4, 1275-1282. Information System and Analytic Hierarchy Wilson, J. P. and Gallant, J. C., 2000, Digital Process. Journal of Cleaner Production, Vol. Terrain Analysis. (New York: Wiley Publisher). 172, 2485-2502. Yeh, H. F., Cheng, Y. S., Lin, H. I. and Lee, C. H., 2016, Mapping Groundwater Recharge Potential Pham, B. T., Jaafari, A., Prakash, I., Singh, S. K., Zone Using a GIS Approach in Hualian River, Quoc, N. K. and Bui, D. T., 2019, Hybrid Taiwan. Sustainable Environment Research, Computational Intellegence Models for Vol. 26, 33-43. Groundwater Potential Mapping. Catena, Vol. Yousif, M., Sabet, H. S., Ghoubachi, S. Y. and Aziz, 182, 104101. A., 2018, Utilizing the Geological Data and Remote Sensing Applications for Investigation Rajasekhar, M., Raju, G. S., Raju, R. S., of Groundwater Occurrences, West El Minia, Ramachandra, M. and Kumar, B. P., 2018, Data Western Desert of Egypt. NRIAG Journal on Comparative Studies of Lineaments Astronomy Geophysics, Vol. 7, 318-333 Extraction from Aster Dem, SRTM, And Cartosat for Jilledubanderu River Basin, Anantapur District, A.P, India by Using Remote Sensing and GIS. Data in Brief, Vol. 20, 1676- 1682. Razandi, Y., Pourghasemi, H. R., Neisani, N. S. and Rahmati, O., 2015, Application of Analytical Hierarchy Process, Frequency Ratio, and Certainty Factor Models for Groundwater International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Wildfire Visualization Time Series from 2014-2019 in Phayao Province, Thailand Jeefoo, P. Geographic Information Science, School of Information and Communication Technology, University of Phayao, 19 Moo 2, Maeka, Muang, Phayao 56000, Thailand E-mail: [email protected], [email protected] Abstract Wildfire prevention comprises of all the activities that can be performed in order to reduce or prevent wildfire occurrences as well as to minimize wildfire potential impacts and damage. The aim of this research was presented wildfire visualize time series from 2014-2019 using time manager plugin in QGIS software. The study area is Phayao province, that highest of wildfire in northern part of Thailand. Time manager plugin was used to visualize time series during the years 2014 to 2019 with 1,993 hotspots especially 2019 year with 582 hotspots. The highest was shown in Pong district with 779 hotspots (39%) of wildfire during March and April of every year. 1. Introduction Free and Open Source Solution for Geoinformatics Wildfire monitoring and management in northern (FOSS4G) provide an alternative approach to part of Thailand is of paramount importance and geospatial software development. FOSS4G has every summer massive forest wildfires break out in become popular as a practical alternative for several areas, leaving behind severe destruction in developers and users and has been successfully used forested and agricultural land, infrastructure and in many applications. The FOSS4G software (QGIS, private property, and losses of human lives. SAGA GIS, GRASS GIS and so on) stack Remotely sensed images have a considerable comprises of system software, data processing tools, potential for mapping burnt areas and fire regimes, data delivery and user interface tools for both particularly at the regional level. This technology is desktop and web-based environment (Choosumrong useful in comparison to other conventional et al., 2016). approaches because of its accuracy, repeatability, and speed in data acquisition, longer historical data Thailand currently faces annual air emission and ease of combination with other thematic data. from forest fires during January to May with impact Satellite-based burnt area mapping is performed on respiratory health and vision, especially in the using different sensors (various spatial and temporal upper northern region of the country (Pollution resolutions) and classification procedures (Daldegan Control Department, 2012 and Junpen et al., 2013). et al., 2014). The main reasons for setting fires include the gathering of forest non-timber products, e.g., honey, A literature review on the surface fire spread mushrooms, bamboo, fuel wood, hunting, models reveals that fire models can develop from agricultural residue burning for land clearing before two types of fuel: grassland and shrub land. The fire the next crop, incendiary fires, and others. Nearly all spread rate depends on three major environmental forest fires in Thailand are classified as surface fires weather parameters, i.e., wind speed at ground level which combust only surface fuel including twigs, (Burrows et al., 2006 and Boboulos and Purvis, dead leaves, dead plants, grass, and undergrowth 2009), moisture control of fuel (Beaza et al., 2002), (Forest Fire Control Division, 2011). All fires are and moisture content of terrain (Fernandes, 2001). human caused, especially by rural people living near Topographic factors include slope of terrain forested areas. Over the period from 1998-2010, the (Boboulos and Purvis, 2009) and fuel total frequency of forest fires was approximately chalacteristics, e.g., high vegetation, vegetation 80,767 events corresponding to a burned surface cover, fuel load in the unit area and fuel size (Bilgili area of 265,472 ha. Approximately 50,872 or 63% and Saglam, 2003) of all forest fire events occurred in the Northern part of Thailand and then over a burned area of Recent advances in Geographic Information approximately 122,688 ha, or 46% of all burned System (GIS) along with Remote Sensing (RS) have areas. Approximately 90% of all forest fire added new dimension to spatial statistics analysis in disaster studies (Jeefoo, 2012). The advances in International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

occurrences took place in deciduous forests The eastern part is mountainous, with an average characterized in to two forest types, mixed height of more than 300 m above sea level. deciduous forests and dry dipterocarp forests (Forest Fire Control Division, 2011). About 28% of the population are predominantly involved in agricultural activities that take place in Remote Sensing will provide objective and an extension of approximately 1,852 square reliable information that could be useful for solving kilometers, including paddy fields, fields crop, the burnt detection (Zheentaev, 2016). Regarding rambutan, and lychee. The average yearly rainfall is fire types of behavior in deciduous forests, the approximately 1,262.40 mm, the highest average Forest Fire Control Division, Thailand, has reported monthly rainfall is September with 232.50 mm, and a surface fire spread rate ranging from 1.70-3.40 lower average monthly rainfall is January with 2.50 m/min which is 1.00-4.00 m of the flame height and mm (Phayao Provincial Agriculture and 110-250 kW/m of fireline intensity. The fire spread Cooperatives Office, 2019). rate might be increased when forest fires occur in steep slope areas. The fire spread rate is a key 2.2 Wildfire Hotspots Data component in planning and decisions for conducting Hotspots burnt detection data, ESRI shape file, prescribed fires and suppressing forest fires. during the years of 2014 to 2019 were collected However, fire spread rate prediction in the area of from NASA via FIRMS website the studies is unknown (Pastor et al., 2003 and [https://firms.modaps.eosdis.nasa.gov]. The data Forest Fire Control Division, 2011). was obtained from the NASA website, with regard to the number of reported wildfire hotspot and Burnt area maps generated from satellite sensors confirmed hotspot locations with geometric with high temporal-resolution and coarse spatial- correction per 24 hour, 48 hour, and 7 day. After the resolution (approximately 1 km.) such as System first wildfire hotspots detection were confirmed, all Pour 1’Observation de la Terra (SPOT- hotspots that had showed in the Phayao province VEGETATION), Advanced Very High Resolution with the following condition: brightness temperature Radiometer (AVHRR), Along Track Scanning ≥ 295 Kelvin, confidence ≥ 60 %, were considered Radiometer (ATSR-2), Geostationary Operational as expects of wildfire hotspot. Data represented only Environmental Satellite (GOES) and Moderate the hotspots detection in the South East Asia region Resolution Imaging Spectrometer (MODIS), fail to and have to process using spatial overlay analysis detect the majority of small and fragmented burnt only Phayao province, Thailand by using QGIS areas (Schroeder et al., 2008, Simon et al., 2004, open source software. Silva et al., 2005 and Daldegan et al., 2014). This study aimed at providing useful information on The data has uploaded by separate comprised of wildfire hotspots detection and time series three categories (Shapefiles, Google Earth KML, management in Phayao province, northern part of Text File (CVS)). A dialog of those data and it can Thailand, using QGIS, Time Manager plugin, free be download by MODIS 1km, VIIR 376m / S-NPP, and open source software and mapping their or VIIRS 375 / NOAA-20 respectively. This study patterns and dynamics of wildfire burning. was selected shapefiles in South East Asia by 24h. The attribute data of hotspots burnt detection was 2. Materials and Methods shown in Table 1, including latitude, longitude, 2.1 Study Area brightness, scan, track, acquisition date, satellite, Phayao province, a province in the northern part of confidence, version, brightness temperature 31 Thailand, had the high burnt area in the northern in (Kelvin), fire radiative power (MW - megawatts) 2019, with 5.16 square kilometers of damaged burnt and day or night. areas, therefore Phayao province was selected as the study area (Figure 1) because of the highest 2.3 Method population ethnics living in the mountain and several types of landscape. Phayao province consists Various software, namely QGIS of 9 districts, 68 sub-districts, and 780 villages. The province is located 733 km north of Bangkok, and (https://qgis.org/en/site/) was used in this study. covered as area of 6,335 square kilometers with geographical location between 18°470N to Meanwhile, Time manager plugin have installed. 19°480N and 100°530E to 100°360E. The province has a population of about 475,215 people Time manager setting, hotspots detection was (NSO, 2018). The western part of the province is the low water resource of Kwan Phayao lake. calculated for each year from 2014 to 2019, the layer name was generated for each year. The ACQ_DATE (Acquisition Date) field have considered for start and end date of the time series (Figure 2). International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Figure 1: Study area, Phayao Province Figure 2: Time manager plugin in QGIS 3. Results again a regular decrease was seen in the number of Hotspot of wildfire detected in most area of Phayao hotspots detection until 2019. The lowest province, causing severe make a living problems occurrence was in 2018 (140 points). As shown in and financial tolls on the population affected. In Table 2, in total 1,993 points were reported. 2019, the total number was 582 points, which is the Wildfire hotspots detection based on month of the highest recorded of hotspot detection for the current hotspots was also determined. The highest hotspot decade. After 2014 (459 points), a gradual decrease was in March with a percentage of 45.01% (897 was seen in the number of hotspots detection until points), while hotspot in the April was 30.96% (617 2019 (in 2016, increase again with 408 points), and points) (Table 3). International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Table 1: FIRM data description Item Detail Description Latitude Longitude Latitude Center of 1km fire pixel but not necessarily the actual Brightness location of the fire as one or more fires can be detected Scan within the 1km pixel. Track Acq_Date Longitude Center of 1km fire pixel but not necessarily the actual Acq_Time location of the fire as one or more fires can be detected Satellite within the 1km pixel. Confidence Brightness temperature Channel 21/22 brightness temperature of the fire pixel Version 21 (Kelvin) measured in Kelvin. Bright_T31 Along Scan pixel size The algorithm produces 1km fire pixels, but MODIS FRP pixels get bigger toward the edge of scan. Scan and Daynight track reflect actual pixel size. Along Track pixel size The algorithm produces 1km fire pixels, but MODIS pixels get bigger toward the edge of scan. Scan and track reflect actual pixel size. Acquisition Date Data of MODIS acquisition. Acquisition Time Time of acquisition/overpass of the satellite (in UTC). Satellite A = Aqua and T = Terra. Confidence (0-100%) This value is based on a collection of intermediate algorithm quantities used in the detection process. It is intended to help users gauge the quality of individual hotspot/fire pixels. Confidence estimates range between 0 and 100% and are assigned one of the three fire classes (low-confidence fire, nominal-confidence fire, or high-confidence fire). Version identifies the collection (e.g. MODIS Collection 6) and source of data processing: Near Real- Time (NRT suffix added to collection) or Standard Version (Collection and Processing (collection only). source) \"6.0NRT\" - Collection 6 NRT processing. \"6.0\" - Collection 6 Standard processing. Find out more on collections and on the differences between FIRMS data sourced from LANCE FIRMS and University of Maryland. Brightness temperature Channel 31 brightness temperature of the fire pixel 31 (Kelvin) measured in Kelvin. Fire Radiative Power Depicts the pixel-integrated fire radiative power in MW (MW - megawatts) (megawatts). Day or Night Day (D), Night (N) International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

2014-03-15 2015-03-15 2016-03-15 2017-03-15 2018-04-15 2019-03-15 Figure 3: Hotspot time series from 2014 - 2019 Table 2: Number of hotspot detection from the year of 2014 – 2019 in Phayao province Year Number of hotspot detection Percent (%) 2014 459 23.03 13.15 2015 262 20.47 7.12 2016 408 7.02 29.20 2017 142 100.00 2018 140 2019 582 Total 1,993 Source: Fire Information for Resource Management System, NASA International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Table 3: Number of hotspot detection by month from 2014 – 2019 Mount Number of Hotspot Detection in Phayao province, Thailand % Total 2014 2015 2016 2017 2018 2019 3.31 66 January 14 8 16 2 12 14 7.88 157 45.01 897 February 54 21 23 11 7 41 30.96 617 3.61 72 March 278 135 135 71 26 252 0.25 0.05 5 April 85 76 211 20 37 188 0.05 1 0.05 1 May 6 5 14 10 3 34 0.00 1 3.11 0 June 2 0 0 1 0 2 5.72 62 114 July 0 0 0 0 0 1 August 000100 September 0 1 0 0 0 0 October 000000 November 1 4 2 0 27 28 December 19 12 7 26 28 22 Source: Fire Information for Resource Management System, NASA The wildfire hotspots detection was illustrated by point for each day during the complete burn, and interpolating the value over the space, using Kernel MODIS sensor was confirmed that hotspots were Density Estimation as interpolation method, as spread all around the province showing a hotspot shown in Figure 3. Maps of monthly hotspots with red colored areas in the west and north of the detection were generated, which indicated the province. This apparent cluster was due to a dynamics of wildfire hotspots distribution through notification effect in the native location, where time during 6 years for the month of March and spatial resolution of hotspots was lower. The highest April (2018) for the years of 2014 to 2019. The brightness temperature was shown 405.7 Kelvin number of hotspot detection by month showed that (date: 2019-04-15) in Pong district and lowest the spatial distribution of wildfire hotspot was brightness temperature was shown 301.2 Kelvin clustered (red color). Maps showed that the hotspot (date: 2019-03-29) in Chiang Kham district. occurred everywhere in the province, even in area in the northern part of the province, which is more In Pong, there was 779 points of hotspot rural with a lower density of villages (Figure 3). detection, affecting 39.09% of the total summary hotspot. Approximately 58.01% (Chiang Kham and Time manager for each year represented the Pong districts) of the occurred in the west of absolute global position and created brightness province, with the highest number occurring in Pong temperature trend in direction and extent for all district and the second highest in Chiang Kham wildfire radius area. Brightness temperature of each district (377 points). While Thu Kham Yaow had location was observed mainly in agriculture areas the lowest number with only 40 points (2.01%). The (Mae Jai and Pong districts) of the Phayao province. summary of hotspot from 2014 – 2019 by district is The purpose of the map was to compare the global shown in Table 4 and Figure 6. distribution of wildfire brightness temperature for 6 serial years. The global pattern was observed almost 4. Discussion same in each year. Maps of monthly hotspots were When it comes to wildfire management, GIS generated, which indicated the dynamics of wildfire provides accurate and comprehensive information to brightness temperature through time during January analyse, visualize, and prioritize values at risk, such to July for the year 2019. A total of 532 hotspots as housing developments, utility infrastructure, detection and wildfire brightness temperature were wildfire, and natural or cultural resources. Within recorded and geo-referenced (Figure 4). The highest the scope of wildfire management, wildfire number of hotspots per day were 32 hotspots (date prevention comprises of all the activities that can be 2019-03-12) in Chiang Kham district. performed in order to reduce or prevent wildfire occurrences as well as to minimize wildfire Figure 5 represents the results for 2019 by potential impacts and damage. For the purpose of month. The hotspots in Phayao province spread mitigation, a GIS plays a vital role in mapping fire- rapidly in all the study area during January and the prone areas, monitoring fuel load and risk wide spatial hotspot distribution was conserved modelling; for preparedness and response it helps in during the peak of the epidemic, at February to fire detection, predicting spread/direction of fire, April 2019. A hotspots detection map was built early warning and coordinate fire-fighting efforts. from the brightness temperature (Kelvin) of each International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Figure 4: Wildfire time series in January 2019 (532 points) January 2019 February 2019 March 2019 April 2019 Figure 5: Wildfire time series in the year of 2019 with totally 532 of hotspots (cont. next page) International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

May 2019 June 2019 July 2019 Figure 5: Wildfire time series in the year of 2019 with totally 532 of hotspots For recovery a GIS can be used as a damage April was 31.76% respectively. In Thailand, assessment tool useful for mapping the extent of wildfire is occurred in summer season especially burn, understanding biological responses due to February to May. differential surface heating (i.e. fire severity) and quantifying extent and pattern of burned areas The wildfire hotspots detection by using (SANPA, 2012). This article aimed at providing brightness temperature (Kelvin) also plays useful information on wildfire hotspots detection important role and it was seen that wildfire hotspots and time series management using QGIS and is generally prevalent in the Phayao province during mapping their patterns and dynamics of wildfire the months of March to April 2019. The highest burning. Time manager plugin proved to be a brightness temperature was shown 405.7 Kelvin valuable tool to analyze the space-time change over (date: 2019-04-15) in Pong district and lowest time. brightness temperature was shown 301.2 Kelvin (date: 2019-03-29) in Chiang Kham district. The study revealed useful information on Approximately 323.20 Kelvin is averaged of brightness temperature time series during 2014 – brightness temperature in the province during 2014 2019 risk assessment to wildfire. In 2019, the total – 2019. However, the limitation in the study was the number was 532 points, which is the highest wildfire hotspots detection data. Due to boundary of recorded of hotspot detection for the current decade administrative line, some hotspots point near the and the lower in 2018 with hotspots detection 140 territorial of Phayao province were moved and points. Moreover, the highest hotspot was in March wildfire hotspots data for these areas was not with a percentage of 46.17%, while hotspot in the available for small area. International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Table 4: Summary of hotspot from 2014 – 2019 by district District Summary of hotspot from Percent (%) 2014 - 2019 Chiang Kham 377 18.92 Chiang Muan 217 10.89 Dok Kham Tai 112 5.62 Jun 96 4.82 Mae Jai 46 2.31 Muang 269 13.50 Pong 779 39.09 Thu Kham Yaow 40 2.01 Thu Sang 57 2.86 1,993 100.00 Total Figure 6: Wildfire hotspot by district from 2014 – 2019 5. Conclusions each year very important to predict burn area. The The results showed that proposed methods and tools methodology is based on notions on general can be beneficial for forest fire control officers to principles of geostatistics and has the potential for visualize and appreciate the distribution and trends application for other disaster. of diffusion shapes of hotspots and to prepare warning and awareness to the atmosphere. Wildfire Acknowledgements hotspots detection diffusion patterns by using I would like to express my sincere gratitude to brightness temperate (Kelvin) may provide useful School of Information and Communication information to support forest fire control officers to Technology (ICT), University of Phayao, Thailand control and predict wildfire spread over hotspot for the support given and providing financial areas only rather than for a whole province. Forest support to this research. fire control officers may service the model to plan an approach to control hotspot by the information References received on supply and hotspots for various day or week. This may save cost and time and make Beaza, M. J., Luis, M. De., Raventos, J. and Escarre, Ministry of Natural Resources and Environment A., 2002, Factors Influencing Fire Behavior in efforts more efficient. In future it would be Shrublands of Different Stand Ages and the important to have regular daily analysis for several Implications for using Prescribed Burning to years to converge faster at hotspots locations and be Reduce Wildfire Risk. J. Environ. Manage. Vol. equipped for deterrent trials. The criteria were 65(2), 199-208. acquisition date and time of wildfire detection of International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Bilgili, E. and Saglam, B., 2003, Fire Behavior in Junpen, A., Garivait, S., Bonnet, S. and Maquis Fuels in Turkey. Forest Ecology and Pongpullponsak, A., 2013, Fire spread prediction Management, Vol. 184(1-3), 201-207. for deciduous forest fires in Northern Thailand. ScienceAsia, Vol. 39, 535-545. Boboulos, M. and Purvis, M. R. I., 2009, Wind and Slope Effects on ROS during the Fire NSO (National Statistical Office), Department of Propagation in East Mediterranean Pine Forest Provincial Administration, Thailand, 2018. Litter. Fire Saf. J. Vol. 44, 764-769. Pastor, E., Zarate, L., Planas, E. and Arnaldos, J., Burrows, N. D., Ward, B., Robinson, A. D. and 2003, Mathematical Models and Calculations Behn, G., 2006, Fuel Dynamics and Fire Systems for the Study of Wildland Fire Behavior in Spinifex Grasslands of the Western Behavior. Prog. Energy Combust. Sci., Vol. 29, Desert. In: Bushfire Conference 2006, Brisbane, 139-153. Australia. Pollution Control Department, 2012, Thailand State Choosumrong, S., Raghavan, V., Jeefoo, P. and of Pollution Report 2010. Pollution Control Vaddadi, N., 2016, Development of Service Department Bangkok, Thailand. [Available Oriented Web-GIS Platform for Monitoring and online at infofile.cd.go.th/mgt/Report_Eng25- Evaluation Using FOSS4G. International 53.pdf] Journal of Geoinformatics, Vol. 12(3), 67-77. Schroeder, W., Prins, E., Giglio, L., Csiszar, I., Fernandes, P.A.M., 2001, Fire Spread Prediction in Schmidt, C., Morisette, J. and Morton, D., 2008, Shrub Fuels in Portugal. Forest Ecology Validation of GOES and MODIS Active Fire Management. Vol. 144, 67-74. Detection Products Using ASTER and ETM+ Data. Remote Sens. Environ., Vol. 112, 2711- Forest Fire Control Division, 2011, Forest Fire 2726. Description. National Park, Wildfire and Plant Conservation Department, Thailand. [Available Silva, J.M.N., Sá, A.C.L. and Pereira, J.M.C., 2005, online at www.dnp.go.th/forestfire/2546/firestat- Comparison of burned area estimates derived istic%20Th.htm, in Thai]. from SPOT-VEGETATION and Landsat ETM+ Data in Africa: Influence of Spatial Pattern and Daldegan, G. A., de Carvalho Junior, O. A., Vegetation Type. Remote Sens. Environ. 96, Guimaraes, R. F., Gomes, R. A. T., Ribeiro, F. 188-201. F. and McManus, C., 2014, Spatial Patterns of Fire Recurrence Using Remote Sensing and GIS Simon, M., Plummer, F., Fierens, F., Hoelzmann, J., in the Brazilian Savanna: Serra do Tombador Arino, O., 2004, Burnt Area Detection at the Nature Reserve, Brazil. Remote Sensing, Vol. 6, Global Scale Using ATSR-2: The GLOBESCAR 9873-9894. Products and Their Qualification. J. Geophys. Res., Vol. 109, 1-16. Jeefoo, P., 2012. Space-Time Analysis Tools of Dengue Epidemics in Chachoengsao province, South African National Space Agency (SANPA), Thailand. International Journal of 2012, Remote Sensing in Disaster Management, Geoinformatics, Vol. 8(3), 9-13. African Leadership Conference. Zheentaev, E., 2016, Application of Remote Sensing Technologies for the Environmental Impact Analysis in Kumtor Gold Mining Company. International Journal of Geoinformatics, Vol. 12, 4, 31-39. International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Aboveground Biomass Estimation in an Upland Tropical Evergreen Forest Environment from Landsat 7 ETM+ Suwanprasit, C.1* and Strobl, J.2 1Department of Geography, Faculty of Social Sciences, Chiang Mai University, 239 Huay Keaw Road Mueang, Chiang Mai, 50200, Thailand, E-mail: [email protected] 2Interfaculty Department of Geoinformatics – Z_GIS, University of Salzburg, Hellbrunnerstrasse 34, Salzburg 5020, Austria, E-mail: [email protected] *Correspondance Author Abstract The study area is an evergreen forest in Phu Hin Rong Kla National Park, Thailand.Tropical evergreen forests are important as they possess the highest biomass amongst tropical forests. Remote sensing has been used in many studies to estimate the volume of aboveground biomass. The gap filled Landsat 7 ETM+ data and the forest observed parameters have been analyzed to find out the quantity of aboveground biomass and carbon sequestration. The forest parameters were measured in 30 randomly selected sample plots measuring 30 m x 30 m. The measurements were then used for computing the aboveground biomass using two allometric equations (TM51 and ND71). The analysis showed the best model for aboveground biomass estimation was a combination of TM51 and ND71 (R2 = 0.658) and the total of aboveground biomass and carbon sequestration were 112,062,010 ton and 5,603,100 ton carbon respectively. The application of the study would be useful for understanding the terrestrial carbon dynamics and global climate change. 1. Introduction estimating biomass including, tree height, basal Estimation of Aboveground Biomass (AGB) is an area, and stand structure, bole diameter at breast essential aspect of studies of carbon stocks and the height (DBH) (Brown, 1997). A strong relation of effects of deforestation and carbon sequestration on biomass with basal area has been found in many the global carbon balance (Alban et al., 1978). It is studies (Rai and Proctor, 1986). There are 3 main also valuable information for many other global approaches to estimating AGB based on (1) Field issues such as forecasting ecosystem productivity, measurement, (2) GIS-based and (3) Remote carbon budget, carbon cycles, nutrient allocation, sensing data (Lu, 2006). and fuel accumulation in terrestrial ecosystems (Amichev et al., 2008, Zhao et al., 2012, Moser et Remote sensing is increasingly used as a source al., 2011, Stape et al., 2008 and Zhang et al., 2007). of data for the efficient and sustainable use of Moreover, for tropical forests, regional change in natural resources and it is considered to be the most biomass is associated with important components of reliable means for spatial biomass estimation for climate change (Brown, 1997). Spatial estimates of tropical regions (Roy and Ravan, 1996). Remote the biomass and volume of trees are becoming more sensing data is often used as the dependent variable and more important to forest management and and therefore there is need to have ground based planning. Unfortunately, studies in tropical forests field surveys as the independent variable (Hahn, are more demanding and time-consuming than in 1984 and Smith, 1986). A number of studies have temperate forests because of high plant species demonstrated the potential of remote sensing for diversity, poorly known plant taxonomy and many biomass studies (Gao et al., 2013, Main-Knorn et of the remaining tropical forests are in remote areas al., 2011, Tian et al., 2012, Miettinen and Liew, with difficult access. Therefore, it is rarely 2009 and Xu et al., 2010). Past studies have shown practicable to collect sufficient field data that covers varying degrees of success in estimating forest the area of interest. Hence, interpolation is often biomass from remote sensing data in temperate and made between widely scattered field sampled tropical forests worldwide (Drake et al., 2003, Benie points. In order to estimate the AGB, forest stand et al., 2005, Anaya et al., 2009, Dubayah et al., parameters such as tree height, basal area, stand 2010, Gleason and Im, 2011, Zheng et al., 2012, structure, bole diameter at breast height (DBH) are Riegel et al., 2013). Remote sensing data is required as independent variables for the models. becoming the primary source for biomass estimation There are a number of commonly used variables for (Lu, 2006). Research on remote sensing based International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

biomass estimation approaches and discussion of Malaysia and Thailand using regression and neural existing issues influencing biomass estimation are networks and discussed the transferability of the valuable for further improving biomass estimation estimates among the regions. Thenkabail et al., performance. (2004) compared narrowband hyperspectral Hypersion data with broadband hyperspatial The potential of remote sensing for estimating IKONOS data, multispectral Advanced Land Imager biomass was reviewed in (Lu, 2006 and Koch, (ALI) and Landsat 7 ETM+ data by using 2010). New remote sensing technology such as classification of complex rainforest vegetation in LiDAR (light detection and ranging), polarimetric southern Cameroon. radar interferometry and hyperspectral data can better estimate forest biomass than optical remote Thus, there is great interest in estimating the sensed data however methodologies and operational biomass of forests and their role in regulating the cost should be concerned (Koch, 2010). Therefore cycling of carbon and nutrients. The objective of optical sensor data would be an alternative source this study was to estimate the AGB and carbon which can directly estimate using several sequestration in an Indo-Malay tropical evergreen approaches such as multiple regression analysis, K forest using field measured parameters and nearest - neighbor, and neural network and vegetation indices derived from Landsat 7 ETM+ indirectly estimated from canopy parameters using Gap-filled spectral reflectance of Phukin Rong Kla multiple regression analysis or various canopy National Park in Central Thailand. reflectance models. Vegetation indices have been successfully used for aboveground biomass 2. Study Area estimation (Lu, 2006). The study area was tropical evergreen forest in Phu Hin Rong Kla National Park. The park is located in Overall, spectral indices of vegetation derived northeast part of Thailand, in Phitsanulok and Loei from the near-infrared and visible (usually red) Provinces. It lies between 16° 53 to 17° 07 E and bands of the remote sensing data are widely 101° 56 to 101° 06 N (Figure 1) and covers an area employed as measures of green vegetation density of 307 km2. The general topography of the park is (Steven et al., 2003). Many vegetation indices have steeply mountainous. The northern part of the park been developed and applied to biophysical in Chaiburi district borders Laos while the southern parameter studies but not all vegetation indices are part runs into Phetchabun Province. The mountain significantly correlated with aboveground biomass range includes the peaks of Phu Phangma, Phu (Lu, 2006). For example, Sader et al., (1989) found Lomlo, Phu Hin Rong Kla and Phu Man Khao that Normalized Difference Vegetation Index which is the tallest peak with an altitude of 1,820 m (NDVI) is not a good predictor of stand structure above sea level. The second tallest is Phu Lomlo variables (e.g. height, diameter of main stem) or with an altitude of 1,664 m. The park is the source total biomass in uneven age, mixed broadleaf forest. of many streams, including Huai Mueat Don, and Despite the application of vegetation indices for Huai Luang Yai (Division, 2002). biomass studies, the accuracy of these biomass models has been questioned because in reality, Phu Hin Rong Kla’s climate is cool the year optical remote sensing provides information on around, especially in the cool season, when canopy leaf density rather than on biomass (Zhang temperature can occasionally drop to freezing point. and Fu, 1999). Mist cover is frequent and the maximum temperature is below 25 oC (Division, 2002). The Several studies have already been carried out to park is covered by mixed deciduous, dry evergreen predict forest attributes using remote sensing and hill evergreen forests. The most important tree imagery. Sader et al., (1989) assessed the feasibility species found in the park are Lagerstroemia of detecting tropical forest successional age class calyculata, Pterocarpus macrocarpus, Shorea and total biomass differences using Landsat TM in siamensis, and S. obtuse, Afzelia xylocarpa the mountain forest of Puerto Rico. National (Division, 2002). Oceanic and Atmospheric Administration (NOAA) Advanced Very High Resolution Radiometer 3. Methods (AVHRR) data for estimating of tropical forest In this study, forest parameters namely diameter at biophysical properties in south - west Ghana was breast height (DBH), tree height, species and used in (Foody et al., 2003). Steininger (2000) also number of trees in each sample plot were measured. tested estimation of aboveground biomass of These parameters were used to estimate the AGB. tropical secondary forest from canopy spectral Remotely sensed data on the other hand was used to reflectance using satellite optical data. (Foody et al., provide vegetation indices for correlation analysis 2003) compared the estimation of tropical forest with the estimated AGB. biomass from Landsat TM data for sites in Brazil, International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Figure 1: Location of study area, Phu Hin Rong Kla National Park, Thailand 3.1 Field Data Collection and Allometric Equations Where, Ts_TB is the Tsutsumi allometric equation, In April 2009, 30 sample plots were randomly Ws is the stem biomass (kg), Wb is the branch selected and forest parameters were measured biomass (kg), Wl is the leaf biomass (kg) and H is (Thompson, 1990). Each plot was composed of an the tree height (cm). area of 30 m x 30 m which is the same as the Landsat 7 Enhanced Thematic Mapper Plus (ETM+) 3.2 Data and Image Processing pixel size. All trees more than 4.5 cm in DBH Two Landsat ETM+ images, covering the study (Viriyabuncha et al., 2002) were measured. The areas and located on path 129 and row 48 were used for analysis. These Landsat scenes were acquired on locations of each plot were determined by the global February 09, 2002 and February 12, 2009. Band 1-5 positioning system (GPS). The measured parameters and 7 of both Landsat 7 ETM+ were used in this were used for estimating the AGB based on study. The 30 m pixel resolution images were geo- allometric relationships. The biomass of individual referenced to UTM coordinates using existing trees was estimated using two regression equations Landsat 7 ETM+ gap-filled image, atmospheric and 1 and 2-5 developed by (Brown,1997) and topographic corrections (by ASTER GDEM 30m (Tsutsumi et al., 1983) respectively. resolution) were also made. Missing data resulting from the failure of the Scan Line Corrector (SLC) Br_TB = 42.69 - 12.800D + 1.242D2 on the Landsat 7 ETM+ sensor, was filled using the Equation 1 Landsat Gap Fill function (Scaramuzza et al., 2004). The boundary of Phu Hin Rong Kla National Park Where, Br_TB is the Brown allometric equation and was then overlaid to exclude pixels outside the study D is the DBH (cm), area. Ws = 0.0509D2H0.919 Equation 2 The radiance and reflectance of Landsat 7 ETM+ Wb = 0.00893D2H0.977 Equation 3 bands were used for calculating the RVI, NDVI, SAVI, ARVI, MSAVI, ND41, ND51, ND53, TM41, Wl = 0.0140D2H0.669 TM51, TM53, VIS123 and MVIS vegetation indices and band ratios (Table 1). These vegetation indices Eqaution 4 and band ratios have been shown to be indicators of overall green biomass, canopy closure, tree density, Ts_TB = Ws + Wb + Wl and tree species diversity (Kaufman and Tanré, Equation 5 1992, Schultz and Halpert, 1995, Lawrence and Ripple, 1998 and Gong et al., 2003). International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Table 1: Vegetation indices and band ratios used in this study Vegetation Index Equation Source ������������������ (Jordan, 1969) Ratio Vegetation Index (Rouse et al., (RVI) ������������������ = ������������������ Normalized Difference ������������������ − ������������������ 1973) Vegetation Index (NDVI) (Huete, 1988) ������������������������ = ������������������ + ������������������ Soil-Adjusted Vegetation (������������������ − ������)(1 + ������) (Qi et al., 1994) Index (SAVI) ������������������������ = (������������������ + ������ + ������) (Kaufman and Modified Soil Adjusted Tanré, 1992) Vegetation Index (MSAVI) ������������������������������ = 2������������������ + 1 − √(2������������������ + 1)2 − 8(������������������ − ������������������) (������������������ − 2������ + ������) Atmospherically Resistant Vegetation Index (ARVI) ������������������������ = (������������������ + 2������ − ������) VIS123 ������������������123 = ������1 + ������2 + ������3 (Lu et al., 2004) MVIS ������������������57 ������������������������ = ������������������123 ND41 (������4 − ������1) ������������41 = (������4 + ������1) ND51 (������5 − ������1) ������������51 = (������5 + ������1) ND53 (������5 − ������3) ������������53 = (������5 + ������3) ND71 (������7 − ������1) This study ������������71 = (������7 + ������1) TM41 ������4 ������������41 = ������1 TM51 ������5 ������������51 = ������1 TM53 ������5 ������������53 = ������3 TM71 ������7 ������������71 = ������1 Remark: Where, VIS123 is the linear combination of Landsat 7 ETM+ Band 1, 2 and 3 MVIS is the ratio of Landsat 7 ETM+ between linear transform of middle infrared region bands (Band 5 and 7) and VIS123 NDab (ND41, ND51 and ND71) is the normalized ratio between Band a and b (where, a and b are the number of Landsat 7 ETM+ bands) TMab (TM41, TM51, TM53 and TM71) is the simple ratio between Band a and b (where, a and b are the number of Landsat 7 ETM+ bands) The best vegetation indices were used in this study Maximum Likelihood (ML), Artificial Neural to estimate the amount of AGB in the whole study Networks (ANNs) and Support Vector Machine area by using regression analysis. A correlation with (SVM). The best accuracy to classify evergreen p < 0.05 is taken as significant. Landsat 7 ETM+ for forest was found using SVM in agreement with the study area was classified into 2 groups many publications (Brown et al., 1999, Heikkinen et (evergreen forest and non-evergreen forest) by al., 2011 and Koetz et al., 2008). comparing three image classifiers that included International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

4. Results cm to 137.77 cm respectively, which represented a 4.1 Forest Stand Parameters full range of stand structures that occurred in the Within the 30 randomly selected sample plots study area (Table 2). The average tree height and whose locations were determined using a Handheld DBH were 13.71 m (standard deviation (S.D.) = GPS, a total of 2,200 trees were measured. The 7.84 m) and 18.90 cm (standard deviation (S.D.) = major tree species identified were Lithocarpus 15.48 cm) respectively. The biomass was estimated calthiformis Rehd. Et Wils. and Syzygium cumini using allometric equations from Brown (1997) and Druce. The highest tree and biggest girth measured Tsutsumi et al. (1983). The averages of AGB from was for a Lithocarpus calathiformis Rehd. et Wils., Brown and Tsutsumi allometric equations per 30 x with a height of 50 m and a girth of 4.33 m. The tree 30 m plot were 567.53 kg and 996.97 kg height and DBH ranged from 3 m to 50 m and 3.39 respectively. Table 2: Descriptive statistics of the validation samples for DBH, tree height and biomass Mean Statistics Max Min 137.77 Forest parameters DBH (cm) 18.90 3.39 50.00 Tree Height (m) 13.71 3.00 21,854.00 9.71 58,282.40 AGB (kg) Br_TB* 567.53 7.20 Ts_TB** 996.97 Remark: * Brown and ** Tsutsumi allometric equations Table 3: Pearson’s correlation coefficients for tree height, average DBH and AGB volume using Brown and Tsutsumi allometric equations Forest Parameters AGB Volume Average DBH Tree Height Br_TB Ts_TB Average DBH 1.000 0.711** 0.888** 0.807** Tree Height 0.711** 1.000 0.564** 0.594** Br_TB 0.888** 0.564** 1.000 0.957** Ts_TB 0.807** 0.594** 0.957** 1.000 Remark: ** Correlation is significant at the 0.01 level (2-tailed) (a) (b) Figure 2: Scatter plot describing the correlation between AGB and average with (a) Brown allometric equation and (b) Tsutsumi allometric equation International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

4.2 Correlation between AGB and Forest Stand composition (Lu et al., 2004). The highest Parameters correlation coefficient is found using Band 1 (- In this study, AGB is defined as the biomass of live 0.715) for the Tsutsumi allometric equation. Table 4 trees greater than 4.5 cm in DBH. The AGB volume summarizes the correlation coefficients between the calculated from both equations has a good selected AGB stand volumes and Landsat 7 ETM+ relationship with forest stand parameters (Table 3 spectral reflectance in the study area using both the and Figure 2). In addition, the correlation coefficient Tsutsumi and Brown allometric equations. (r) between average DBH and AGB from Brown allometric equation (r = 0.888) is higher than AGB 4.4 AGB and Vegetation Indices and Band Ratios calculated using the Tsutsumi allometric equation (r Not all vegetation indices are significantly (p < = 0.807) and tree height also has a high correlation 0.05) related to forest stand parameters (Lu et al., with average DBH. The correlation coefficient of 2004). Table 5 and Table 6 summarize the the Brown allometric equation (r = 0.564) is similar correlation coefficients between AGB volume to the correlation derived from the Tsutsumi estimated from Brown and Tsutsumi allometric allometric equation (r = 0.594), even though the equations with vegetation indices in the study area. Brown allometric equation has no tree height The results show that all indices were significantly parameter included in it. correlated with AGB in the sample plots using both equations except MSAVI. The highest correlation 4.3 Aboveground Biomass and Landsat 7 ETM+ coefficient was found using the TM51 with the Spectral Reflectance Tsutsumi allometric equation. Moreover, the All selected AGB stand volumes from Tsutsumi and correlation coefficients for AGB volume that were Brown allometric equations have significant (p < calculated from the Brown allometric equation are 0.05) but negative correlations with Landsat 7 higher than AGB values that were calculated from ETM+ reflectance. Moreover, Band 1 to 3 are the Tsutsumi allometric equation for all indices significantly correlated with AGB using both except RVI. However, most vegetation indices are equations but such relationships vary depending not better related to forest stand parameters than upon the characteristics of the study areas such as ETM+ spectral reflectance except TM51 (Table 6). forest stand density, vegetation age and species Table 4: Pearson’s correlation coefficient from two AGB volumes by using Brown equation and Tsutsumi allometric equations and six Landsat 7 ETM+ spectral bands AGB Landsat 7 ETM+ spectral bands Band1 Band2 Band3 Band4 Band5 Band7 -0.226 ns BR_TB -0.638** -0.589** -0.593** -0.114ns -0.213 ns -0.186 ns TS_TB -0.715** -0.673** -0.606** -0.185 ns -0.188 ns Remark: ** Correlation is significant at the 0.01 level (2-tailed). ns is not significant. Table 5: Pearson’s correlation coefficient from two AGB volumes by using Brown and Tsutsumi allometric equations and vegetation indices Index RVI NDVI SAVI ARVI MSAVI VIS123 MVIS 0.635** Br_TB 0.527** 0.569** 0.569** 0.626** 0.11 ns -0.645** 0.710** Ts_TB 0.496** 0.572** 0.572** 0.646** 0.047 ns -0.706** Remark: ** Correlation is significant at the 0.01 level (2-tailed). ns is not significant. Table 6: Pearson’s correlation coefficient from two AGB volumes by using Brown and Tsutsumi allometric equations and band ratios. Index ND41 ND51 ND53 ND71 TM41 TM51 TM71 TM53 Br_TB 0.659** 0.522** 0.630** Ts_TB 0.627** 0.517** 0.628** 0.399* 0.680** 0.721** 0.599** 0.620* Remark: 0.698** 0.603** 0.663** 0.487** 0.703** ** Correlation is significant at the 0.01 level (2-tailed) * Correlation is significant at the 0.05 level (2-tailed) International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

4.5 Biomass and Carbon Sequestration Estimation carbon/whole area (307 km2). The average of carbon 4.5.1 Model selection storage of hill evergreen forest was 497 ton The top five models with the highest coefficient carbon/ha. Figure 3(b) and Table 8 showed the determination are showed in Table 7. Cubic model carbon sequestration in evergreen forest area in the which added TM51 in the model was the highest study area. coefficient of determination (R2 = 0.501) among models using Brown allometric equation. For the 5. Discussion Tsutsumi allometric equation, the highest coefficient The study found an inverse relationship of spectral of determination was a multiple regression model reflectance with AGB information except Band 4 which included TM51 and ND71 in the model (R2 = where almost no relationship was apparent. Similar 0.658). conclusions were drawn in many of the previous studies (Horler et al., 1983, Spanner et al., 1990, 4.5.2 Aboveground Biomass and Carbon Roy and Ravan, 1996). The inverse relationship Sequestration in the Evergreen Forest Area between the biomass and spectral reflectance could The Tsusumi allometric equation had the highest be the result of increased canopy shadowing within correlation coefficient (R2 = 0.658). The model was larger stands and the decreased understory fitted using multiple linear regression (including brightness (soil brightness) due to increased density TM51 and ND71 as variables in the model). The of biomass (Spanner et al., 1990). There was an estimated AGB in evergreen forest is shown in inverse relationship between amount of shadow and Figure 3(a) and Table 8 presents descriptive reflectance in all bands (Ahlcrona, 1988). Results statistic. The total of AGB in the whole study area from the study by (Horler et al., 1983) specified this was 112,062,010 ton. The average of AGB in phenomenon for particular spectral regions, for evergreen forest areas was 872 ton/ha. The result example band 5 or 7 of Landsat TM and stated that was similar to that found by (Viriyabuncha et al., shadowing is a factor at least as important as leaf 2002) at Doi Suthep - Pui National Park, Chiang moisture content influencing the spectral reflectance Mai, evergreen forest mixed deciduous forest where of forests in the shortwave infrared spectral region. AGB was in the range 32-903 ton/ha. The carbon Thus, the higher spectral reflectance of the sample content estimation would be about 50% of the plots with less biomass can be explained partially by amount of total aboveground biomass (Dixon et al., a lower amount of shadow, which results in a higher 1996). Therefore, the carbon sequestration of contribution to the spectral radiance from the evergreen forest was calculated and the total carbon background soil. in the whole area was stored as 5,603,100 ton Table 7: Top five of the highest coefficient determination between two allometric equations and spectral reflectance Models R2 Brown allometric equation 0.501 1 AGB = 4.024 – 0.010(TM51) + 0.0000226(TM51)2– 0.000000008(TM51)3 0.499 2 AGB = 3.151 – 0.004(TM51) + 0.00000958(TM51)2 0.480 3 AGB = 7.822 – 0.017(TM41) + 0.000056(TM41)2 – 0.000000027(TM41)3 0.477 0.468 4 AGB = 4.904 + 0.003(TM41) + 0.0000126(TM41) 2 0.658 5 AGB = 843.853e-0.002(B1) 0.602 0.573 Tsutsumi allometric equation 0.547 1 AGB = 421.905(TM51) – 1,823.46(ND71) – 204.4 0.546 2 AGB = 126.233(TM51) – 0.364(VIS123) + 880.342 3 AGB = 730.67e-0.087(B1) 4 AGB = 3.124 – 0.002(TM51)+ 0.000004 (TM51)2– 0.000000001(TM51)3 5 AGB = 2,179.93 – 0.7(VIS123) – 0.01(VIS123)2 + 0.000000277(VIS123)3 Table 8: Descriptive statistic of AGB estimated and carbon sequestration from Tsutsumi allometric equation, TM51 and ND71 for study area (307 km2) Descriptive Statistic Minimum Maximum Mean Range Summary (ton/ha) (ton/ha) (ton/ha) (ton/whole area) AGB 34,216 34,216 Carbon sequestration 0 17,108 872 17,108 112,062,010 0 467 5,603,100 International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Figure 3: (a) Aboveground biomass and (b) Carbon sequestration estimated in evergreen forest area Previous studies have shown that there is a high pixels. Moreover (Okuda et al., 2004) concluded correlation between AGB and vegetation indices that Landsat data is probably insufficient to detect (Sader et al., 1989, Huete et a., 1997, Yemefack et local difference in AGB by using only visible, near- al., 2006, Zheng et al., 2004, Anaya et al., 2009 and infrared and shortwave infrared wavelengths (0.4 – Helmer et al., 2009). Conversely, there are some 2.5 mm). In addition, the study area was upland studies that have found low or weak correlation tropical evergreen forest, therefore, there are between AGB and vegetation indices (Foody et al., difficulties in measuring these variables accurately 2003, Lu et al., 2004 and De La Cueva, 2008). In and inter-species differences should ideally be this study, the correlations between some vegetation accounted for throughout. The species composition, indices derived from Landsat 7 ETM+ spectral forest stands structure and associated canopy bands such as Band 4, Band 5 and Band 7 are not shadows, and vegetation vigor, are considered to be significantly correlated with the AGB in the study important factors affecting vegetation reflectance area. The inverse relationship typically (as shown in (Lu et al., 2004). Furthermore, Landsat 7 ETM+ Table 4) can be found between spectral band data data primarily captures canopy information instead and ABG (Roy and Ravan, 1996). They proposed of individual tree information due to its limited that the underlying causes were; (i) increased spatial resolution. For these reasons, future research canopy shadowing within larger stands (ii) may focus on the integration of multi - source data, decreased soil brightness due to increased density particularly RADAR or LIDAR with higher with which biomass increases. penetrating power (Saatchi et al., 2011, Solberg et al., 2014, Montesano et al., 2015), which involves The low correlation was attributed to the the effective integration of remote sensing, GIS and inability of limited data set with such spectral modeling techniques. resolution to account for the variability of forest biophysical features that relate to biomass. 6. Conclusion Therefore, the disadvantage of Landsat data is a This study attempted to provide specifically adapted limitation of resolution. Canopy discontinuities are methods to extract aboveground biomass/carbon sometimes lost in the 30 m x 30 m resolution information from optical remote sensed data from (Feldpausch et al., 2006) because of mixing of International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Indo-Malay tropical forest. The multiple regression Ecology and Management, Vol. 256(11), 1949- models were developed based on integration of 1959. doi: DOI 10.1016/j.foreco.2008.07.020. Landsat 7 ETM+ images and forest stand Anaya, J. A., Chuvieco, E. and Palacios-Orueta, A., parameters. The Tsutsumi allometric equation was 2009, Aboveground Biomass Assessment in found to be the best predictor of AGB from Landsat Colombia: A Remote Sensing Approach. Forest 7 ETM+ data using TM51 and ND71 (R2 = 0.658) Ecology and Management, Vol. 257(4), 1237- in this study. Usually, the Band 1 (Blue region) 1246. doi:DOI10.1016/j.foreco.2008.11.016. shows a very poor relationship between forest Benie, G. B., Kabore, S. S., Goita, K. and Courel, variables and spectral reflectance satellite image M. F., 2005, Remote Sensing-Based Spatio- data. However, the relationship can be significantly Temporal Modeling to Predict Biomass in improved using a correlation of Band 1 with Band Sahelian Grazing Ecosystem. Ecological from the infrared region (Band 4, 5 and 7) to Modelling, Vol. 184(2-4), 341-354. doi: DOI formulate a band ratio that is better correlated with 10.1016/j.ecolmodel.2004.10.012. AGB. Moreover, this study also found that NDVI Brown, S., 1997, Estimating Biomass and Biomass was not the best vegetation index for estimating Change of Tropical Forests: A Primer. FAO AGB in this area. Forestry Paper-134. Rom, Italy: Food and Agriculture Organization of the United Nations. Although optical remote sensed data still Brown, S. L., Schroeder, P. and Kern, J. S., 1999, produce lower accuracy for AGB estimation than Spatial Distribution of Biomass in Forests of the active remote sensed data such as LiDAR and Eastern USA. Forest Ecology and Management, Radar, the data sources availability and longer Vol. 123(1), 81-90. archive are the advantages of optical remote sensed De La Cueva, A. V., 2008, Structural Attributes of data. Therefore AGB estimation methods Three Forest Types in Central Spain and Landsat development in tropical forest using optical images ETM Plus Information Evaluated with is still very important due to the difficulty in Redundancy Analysis. International Journal of gathering ground-truth data representative of a large Remote Sensing. Vol. 29(19), 5657-5676. doi: area. However the results can be used to guide the Doi 10.1080/01431160801891853. selection of suitable Landsat 7 ETM+ band(s) and Division, Forest Land Resources, 2002, Master Plan vegetation indices for AGB estimation in Indo- Database of Phu Hin Rongla National Park Malay tropical forest in the future, the integration of Management, Phitsanulok and Loei Provinces. multisource data for better AGB estimation should In Forest Land Resources Division. Bangkok: be more explored and concerned. Forest Land Resources Division. Dixon, R. K., Krankina, O. N. and Kobak, K. I., 7. Acknowledgments 1996, Global Climate Change Adaptation: Authors wish to gratefully acknowledge support Examples from Russian Boreal Forests. from Dr.Raymond J. Ritchie for his helpful Adapting to Climate Change: An International comments on the paper. This research was Perspective, 359-373. supported by the Austrian agency for international Drake, J. B., Knox, R. G., Dubayah, R. O., Clark, D. mobility and cooperation in education, science and B., Condit, R., Blair, J. B. and Hofton, M., 2003, research (OeAD-GmbH) for field data collecting. Above-Ground Biomass Estimation in Closed Canopy Neotropical Forests Using Lidar Remote References Sensing: Factors Affecting the Generality of Relationships. Global Ecology and Ahlcrona, E., 1988, The Impact of Climate and Man Biogeography, Vol. 12(2), 147-159. doi: DOI on Land Transformation in Central Sudan: 10.1046/j.1466-822X.2003.00010.x. Applications of Remote Sensing, Lund Dubayah, R. O., Sheldon, S. L., Clark, D. B., University Press. Hofton, M. A., Blair, J. B., Hurtt, G. C. and Chazdon, R. L., 2010, Estimation of Tropical Alban, D. H., Perala, D. A. and Schlaegel, B. E., Forest Height and Biomass Dynamics Using 1978, Biomass and Nutrient Distribution in Lidar Remote Sensing at La Selva, Costa Rica. Aspen, Pine, and Spruce Stands on the Same Journal of Geophysical Research- Soil Type in Minnesota. Canadian Journal of Biogeosciences, Vol. 115(G2). DOI: 10.1029- Forest Research, Vol. 8(3). /2009JG000933. Feldpausch, T. R., McDonald, A. J., Passos, C. A. Amichev, B. Y., Burger, J. A. and Rodrigue, J. A., M., Lehmann, J. and Riha, S. J., 2006, Biomass, 2008, Carbon Sequestration by Forests and Soils Harvestable Area, and Forest Structure on Mined Land in the Midwestern and Appalachian Coalfields of the US. Forest International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International

Estimated from Commercial Timber Inventories Huete, A. R., Liu, H. Q. and vanLeeuwen, W. J. D., and Remotely Sensed Imagery in Southern Amazonia. Forest Ecology and Management, 1997.,The Use of Vegetation Indices in Forested Vol. 233(1), 121-132. Foody, G. M., Boyd, D. S. and Cutler, M. E. J., Regions: Issues of Linearity and Saturation. 2003, Predictive Relations of Tropical Forest Biomass from Landsat TM Data and their IGARSS'97. 1997 IEEE International Transferability between Regions. Remote Sensing of Environment. Vol. 85(4), 463-474. Geoscience and Remote Sensing Symposium doi: Doi 10.1016/S0034-4257(03)00039-7. Gao, Y. H., Liu, X. X. Min, C. C., He, H. L., Yu, G. Proceedings. Remote Sensing - A Scientific R., Liu, M., Zhu, X. D. and Wang, Q., 2013, Estimation of the NorthSouth Transect of Vision for Sustainable Development. Vols I- Eastern China Forest Biomass Using Remote Sensing and Forest Inventory Data. International Iv:1966-1968. DOI: 10.1109/IGARSS.1997.60- Journal of Remote Sensing, Vol. 34(15), 5598- 5610. Doi 10.1080/01431161.2013.794985. 9169. Gleason, C. J. and Im, J., 2011. A Review of Remote Sensing of Forest Biomass and Biofuel: Jordan, C. F., 1969, Derivation of Leaf-Area Index Options for Small-Area Applications. Giscience & Remote Sensing, Vol. 48(2):141-170. Doi from Quality of Light on Forest Floor. Ecology, 10.2747/1548-1603.48.2.141. Gong, P., Pu, R., Biging, G. S. and Larrieu, M. R., Vol. 50(4), 663-666. 2003, Estimation of Forest Leaf Area Index Using Vegetation Indices Derived from Kaufman, Y. J. and Tanré, D., 1992, Hyperion Hyperspectral Data. IEEE Transactions on Geoscience and Remote Atmospherically Resistant Previous Sensing, Vol. 41(6), 1355 -1362, DOI: 10.1109/TGRS.2003.812910. Termvegetation Indexnext Term (ARVI) for Hahn, J. T., 1984, Tree Volume and Biomass Equations for the Lake States. Research Paper EOS-MODIS. IEEE Transactions on Geoscience NC-250. St. Paul, MN: U.S. Dept. of Agriculture, Forest Service, North Central Forest and Remote Sensing, Vol. 30(10). DOI: 10.110- Experiment Station. Heikkinen, V., Korpela, I., Tokola, T., Honkavaara, 9/36.134076 E. and Parkkinen, J., 2011, An SVM Classification of Tree Species Radiometric Koch, B., 2010, Status and Future of Laser Signatures Based on the Leica ADS40 Sensor. IEEE Transactions on Geoscience and Remote Scanning, Synthetic Aperture Radar and Sensing, Vol. 49(11), 4539-4551. doi: Doi 10.1109/Tgrs.2011.2141143. Hyperspectral Remote Sensing Data For Forest Helmer, E. H., Lefsky, M. A. and Roberts, D. A., 2009, Biomass Accumulation Rates of Biomass Assessment. ISPRS Journal of Amazonian Secondary Forest and Biomass of Old-Growth Forests from Landsat Time Series Photogrammetry and Remote Sensing, Vol. and the Geoscience Laser Altimeter System. Journal of Applied Remote Sensing, Vol. 3(1). 65(6), 581-590. doi: DOI 10.1016/j.isprsjprs- Doi 10.1117/1.3082116. Horler, D. N. H., Dockray, M. and Barber., J., 1983, .2010.09.001. The Red Edge of Plant Leaf Reflectance. International Journal of Remote Sensing, Vol.4 Koetz, B., Morsdorf, F., van der Linden, S., Curt, (2), 273-288. Huete, A. R., 1988, A Soil-Adjusted Vegetation T. and Allgower, B., 2008, Multi-Source Land Index (Savi). Remote Sensing of Environment, Vol. 25(3), 295-309. Cover Classification for Forest Fire Management Based on Imaging Spectrometry and LiDAR Data. Forest Ecology and Management, Vol. 256(3), 263-271. doi: DOI 10.1016/j.forec- o.2008.04.025. Lawrence, R. L. and Ripple, W. J., 1998, Comparisons among Previous Termvegetation Indicesnext Term and Bandwise Regression in a Highly Disturbed, Heterogeneous Landscape: Mount St. Helens, Washington. Remote Sensing of Environment , Vol. 64(1), 91-102. https://doi.org/10.1016/S0034-4257(97)00171-5. Lu, D. S., 2006, The Potential and Challenge of Remote Sensing-Based Biomass Estimation. International Journal of Remote Sensing, Vol. 27(7), 1297-1328. doi: Doi 10.1080/014311- 60500486732. Lu, D. S., Mausel, P., Brondizio, E. and Moran, E., 2004, Relationships between Forest Stand Parameters and Landsat TM Spectral Responses in the Brazilian Amazon Basin. Forest Ecology and Management, Vol. 198(1-3),149-167. doi: DOI 10.1016/j.foreco.2004.03.048. Main-Knorn, M., Moisen, G. G., Healey, S. P., Keeton, W. S., Freeman, E. A. and Hostert, P., 2011, Evaluating the Remote Sensing and Inventory-Based Estimation of Biomass in the Western Carpathians. Remote Sensing, Vol. 3(7), 1427-1446. doi: Doi 10.3390/Rs3071427. International Journal of Geoinformatics, Volume 16, No. 4, October - December 2020 Online ISSN 2673-0014/ © Geoinformatics International