http://journal.uir.ac.id/index.php/JGEET E-ISSN : 2541-5794 P-ISSN : 2503-216X Journal of Geoscience, Engineering, Environment, and Technology Vol 6 No 4 2021 184 Putra, U.G. et al./ JGEET Vol 6 No 4/2021 RESEARCH ARTICLE Interpretation of Subsurface Fault Through Multi-Level Second Vertical Derivative Gravitational Data in Bittuang Geothermal Working Area, South Sulawesi, Indonesia Ullil Gunadi Putra1*, William Jhanesta2,3, Iskandarsyah1 1 Geophysics Study Program, Department of Geosciences, Universitas Indonesia, Depok 16424, Indonesia 2 Laboratory of Geophysical Modelling, Department of Geosciences, Universitas Indonesia, Depok 16424, Indonesia 3 enerGIS Indonesia, Wisma Anugerah, Depok 16425, Indonesia * Corresponding author : ullil.gunadi@ui.ac.id Tel.:+62-895-6027-23393 Received: Sept 19, 2021; Accepted: Dec 24, 2021. DOI: 10.25299/jgeet.2021.6.4.7744 Abstract The research was conducted in Bittuang, Tana Toraja Regency, South Sulawesi Province, as one of the geothermal prospect areas and targets for the initial stage of the Government exploration drilling program for the 2020-2024 period. One aspect of geothermal is the manifestation control structure as a fluid migration path from below the surface. Therefore, identification of existing struc tures in the Bittuang geothermal area was carried out and confirmed the surface geological structure contained in the Bittuang geothermal geological map. In determining the presence of a fault and knowing its characteristics such as the type of fault, the direction of the d ip, and the magnitude of the dip of the fault, the gravity data is processed using the multi-level second vertical derivative (ML-SVD) method. To strengthen the interpretation, the results from the ML-SVD were matched with the data from the horizontal gradient (HG) method and the geological data of the structure of the study area. From this process, there are 27 faults in the Bittuang geothermal area, two of which are indicated as controlling faults for the manifestation of the Balla group and the Cepeng group. This research is expected to describe faults in the Bittuang geothermal area, which can support detailed exploration activities. Keywords: Gravity, ML-SVD, Structure, Geothermal 1. Introduction Energy needs in Indonesia and the world have increased yearly, while fossil energy is decreasing and is not environmentally friendly, making it an option that new renewable energy (EBT) is needed as an alternative. One of the renewable energy sources that can be utilized which is geothermal energy. Geothermal is a source of thermal energy in hot water, water vapor, rocks, and associated minerals and other gases that are genetically inseparable in a geothermal system. This research was conducted in the Bittuang area, South Sulawesi, with geothermal manifestations in fumaroles, hot springs, alteration of rocks, and dead solfatara fields. The Bittuang field is one of the targets for the initial stage of the exploration drilling program carried out by the Indonesian Government for the 2020-2024 period. In geothermal exploration activities, geophysical methods have a crucial role in determining geothermal aspects, one of which is the gravity method. The gravity method is widely used and applies to detecting subsurface structures associated with geothermal systems (Jhanesta and Supriyanto, 2021a; Nishijima and Naritomi, 2017; Supriyanto et al., 2017; Uwiduhaye et al., 2018). Although gravity is not the main geophysical method used in geophysical exploration, previous studies show that structural analysis with gravity data is quite effective. This study aims to detect the presence of subsurface structures as a controller for the emergence of geothermal manifestations with a physical approach. The ML-SVD method is used to map the distribution of the structure and estimate the dip value and slope of the fault. Finally, this research can know the fault characterization much better so that it is expected to reduce uncertainty in geothermal exploration and increase the probability of success in drilling activities. 2. Field Review Administratively, the Bittuang geothermal area is included in Tana Toraja Regency, South Sulawesi Province. Based on the geological map (Fig 1), the rock in the Bittuang Geothermal area has ten rock units, consisting of seven volcanic rock units, one sedimentary rock unit, one intrusive rock unit, and one metamorphosed rock unit. The stratigraphic order from the youngest age to the oldest age is Lava of Mt Karua 3, Pyroclastic Fall of Mt Karua, Pyroclastic Flow of Mt Karua, Lava of Mt Karua 2, Lava of Mt Karua 1, Intrusion of Rattebombong, Lava of Mt Ruppu, Lava of Mt Panusuk, Sandstone, and Metamorphic rocks (Indonesian Geological Agency, 2009). Meanwhile, the structures found in the Bittuang area are the collapsed structure, normal faults, and several horizontal faults. http://journal.uir.ac.id/index.php/JGEET Putra, U.G. et al./ JGEET Vol 6 No 4/2021 185 Fig 1. Geological map of Bittuang geothermal working area after Indonesian Geological Agency (2009). (BLA = Balla Groups Hot Springs, CPG = Cepeng Groups Hot Springs, FUM = Fumarole) Fig 2. Ternary diagram of Bittuang geothermal manifestations after Kusnadi and Setiawan (2009). It reveals that the Balla groups have chloride characteristics, while the Cepeng groups are bicarbonate types. Collapse structure is a part that has collapsed due to the void in the earth's bowels due to the eruptive activity of Mt Karua. The normal faults found in the Bittuang geothermal area generally have a northeast-southwest, southeast- northwest direction, and some have an almost south-north direction, namely the Tombilangi fault and the Balla fault controls the emergence of the Balla and Cepeng geothermal manifestation. The strike – slip fault in this area has a northeast-southwest direction that cuts rocks and pre- existing structures that cause the shift of rocks around the area (Indonesian Geological Agency, 2009). 186 Putra, U.G. et al./ JGEET Vol 6 No 4/2021 Fig 3. Residual anomaly gravity map. It reveals low anomaly distribution on the center of research area and high anomaly distribution on the south to north-eastern side. Fig 4. (a) Horizontal gradient anomaly map and (b) Second vertical derivative anomaly map of Bittuang geothermal field. According to Kusnadi and Setiawan (2009), the Bittuang Geothermal area has two groups of geothermal manifestations based on the location where these manifestations appear. These groups are the Balla manifestation group and the Cepeng manifestation group. The Balla manifestation group has geothermal manifestations in hot springs (BLA-1 and BLA-2 in Fig 1), fumaroles (FUM in Fig 1), alteration rocks, and inactive or dead solfatara zones. In comparison, the Cepeng manifestation group has geothermal manifestations in the form of hot springs (CPG-1 and CPG-2 in Fig 1). Putra, U.G. et al./ JGEET Vol 6 No 4/2021 187 Fig 5. ML-SVD map of Bittuang geothermal field. The spectrum colours represent the various upward continuation levels. Red solid lines show the slicing lines of HG and ML-SVD. Fig 6. HG and ML-SVD graphs indicate subsurface structure related to the geothermal system. Dash lines represent the HG value, while the solid lines represent the ML-SVD value for each slicing lines. Balla-1, Balla-2, and Balla-3 hot springs which are in the Balla manifestation group have characteristics such as clear water, sulfur smell, slightly salty taste, and bubbling gas, with temperatures between 48.1℃ to 96.7℃ at air temperature 22.5℃ with a pH of about 8.4 and a discharge of 1 liter/second, and around the hot springs there are 188 Putra, U.G. et al./ JGEET Vol 6 No 4/2021 deposits of sintered silica, for more details can be seen in table 1. Table 1 Characteristics of Balla-1, Balla-2, and Balla-3 hot springs after Indonesian Geological Agency (2009). Characteristics Balla-1 Balla-2 Balla-3 Temperature (℃ ) 96.7 90.5 48.1 Air Temperature (℃ ) 22.5 21.6 18.8 pH 8.4 7.75 5.4 Debit (liter/second) 1 0.5 2 Electrical Conductivity (µs/cm) 9700 8400 222 SiO2 (mg/liter) 179.7 153.04 114.85 Na (mg/liter) 1848 15.04 50.80 K (mg/liter) 92.6 68.5 13.61 Ca (mg/liter) 5.57 12.08 9.28 Mg (mg/liter) 1.5 1 0.07 Li (mg/liter) 10.4 7.9 0.04 Cl (mg/liter) 2459.71 2144 80.13 SO4 (mg/liter) 378.58 331.67 4 HCO3 (mg/liter) 223 91.17 49.35 B (mg/liter) 85.18 76.86 0.76 As (mg/liter) 8 331.67 0.1 F (mg/liter) 1 0.5 0.3 Al (mg/liter) 0.16 0.1 - Fe (mg/liter) 0.09 0.07 0.07 CO3 (mg/liter) - - - NH4 (mg/liter) - - - Ion balance (%) 2.43 0.2 -2.07 Meanwhile, the cepeng-1 and cepeng-2 hot springs which are in the cepeng manifestation group have the characteristics such as clear water, sour taste, and no gas bubbles, for more details can be seen in table 2. Table 2 Characteristics of Cepeng-1 and Cepeng-2 hot springs after Indonesian Geological Agency (2009). Characteristics Cepeng-1 Cepeng-2 Temperature (℃) 37.6 39.8 Air Temperature (℃) 22.1 20.1 pH 6.28 5.97 Debit (liter/second) 1 2 Electrical Conductivity (µs/cm) 2630 1241 SiO2 (mg/liter) 159.86 170.25 Na (mg/liter) 210.96 65.68 K (mg/liter) 13.30 4.56 Ca (mg/liter) 190.20 119.34 Mg (mg/liter) 105.60 45.90 Li (mg/liter) 1.20 0.48 Cl (mg/liter) 376.40 105.11 SO4 (mg/liter) 60.33 4 HCO3 (mg/liter) 998.25 592.94 B (mg/liter) 9.21 0.93 As (mg/liter) - - F (mg/liter) 0.5 0.5 Al (mg/liter) - - Fe (mg/liter) 4.73 3.09 CO3 (mg/liter) - - NH4 (mg/liter) - - Ion balance (%) -0.24 -0.58 Based on the chemical characteristics and types of hot water using a triangular diagram of Cl-SO4-HCO3 and Na-K- Mg, presented in Fig 2 (Giggenbach, 1991), Balla-1 hot springs, Balla-2, and Balla-3 are located in the chloride zone. This zone is thought to be related to water originating from the reservoir of the geothermal system. The Balla manifestation group is indicated to rule out near the upflow zone due to the presence of fumaroles, high temperature in hot springs, and silica sinter. While the Cepeng manifestation group is indicated as an outflow zone because it has a low temperature, neutral water pH, bicarbonate water characteristics, and the dominant influence of surface water in the formation of hot water (Kusnadi and Setiawan, 2009). Furthermore, based on the Na-K-Mg ternary diagram, Balla-1 hot springs and Balla-2 hot springs are in the partial equilibrium zone, where this zone indicates that the existing fluid has interacted with the surrounding rock or known as water-rock interaction, thus forming hot water with a high temperature (90-96 ℃). Meanwhile, the Balla-3 hot spring is in the immature water zone, which indicates that the temperature of the manifestation has a lower temperature, which is between 37-48 ℃. The condition is influenced by the fluid that interacts with the surrounding rocks in a hot state and interacts with meteoric water. 3. Methodology This study uses satellite gravity data from the Global Gravity Model plus (GGMplus), which has a resolution accuracy of up to 200 m (Hirt et al., 2013). For gravity processing, GGMplus data is extracted using MATLAB software to obtain the gravity disturbance value. Before carrying out advanced geological interpretation, this data needs to be processed further by performing basic gravity data processing, such as free air correction, Bouguer correction, and terrain correction to obtain a complete Bouguer anomaly (Abdel, 2011). The subsurface delineation process is carried out using derivative analysis, namely horizontal gradient (HG) and ML-SVD. This method is able to detect geological boundaries with a high-value response, but only on structures that are formed vertically (Blakely, 1995). Calculations can be done with the following equation: 𝐻𝐺 = √( ∂𝑔 ∂𝑥 ) 2 + ( ∂𝑔 ∂𝑦 ) 2 (1) Where ∂𝑔 ∂𝑥 and ∂𝑔 ∂𝑦 are the first derivatives of the gravitational anomaly in the x and y directions, respectively. The high anomalous response in HG does not always indicate a structure related to the geothermal system (Jhanesta and Supriyanto, 2021a). Therefore, other methods and supporting data are needed to strengthen the interpretation. The ML-SVD method is a technique that is used to estimate dip and slope faults by applying an SVD filter on the results of upward continuation (Rosid and Naufal, 2019). In this study, the upward continuation process is carried out from 0 m to 2000 m with a continuation interval of 200 m. The dip calculation is done by calculating the arctan value from the linear regression equation gradient between fault location and upward continuation level. At the same time, the slope direction is estimated by analyzing the shifting direction of the zero-SVD value (Jhanesta and Supriyanto, 2021b). 4. Results and Discussion 4.1. Gravity Anomaly Putra, U.G. et al./ JGEET Vol 6 No 4/2021 189 The estimation of the average rock density of the study area using the F-H method, based on the curve generated using this method, the gradient from the line equation is 2.40, which is also the average density value of the study area. Based on geological data, the study area is dominated by pyroclastic flows consisting of ash-lapilli-sized tuff. There are also pumice, sticky, and rhyolite-dacitic compositions. The results of the petrographic analysis conducted by the Indonesian Geological Agency; the rock is rhyolite. Pyroclastic flows are extrusive igneous rocks. According to Schon (2015), the density of extrusive igneous rocks with a rhyolite composition is between 2.15 gr/cc to 2.6 gr/cc. Thus, the average density of the study area of 2.40 gr/cc obtained represents the pyroclastic flow with the dominant rhyolite composition in the study area. The gravity anomaly residual map (Fig 3) shows the high anomaly in red while the low anomaly is light purple with values ranging from -8.5 mGal to 5.7 mGal. If referring to the geological map of the research area (Fig 1), the high anomaly values are in Pleistocene rhyolite rock, which is a pyroclastic flow unit of Mt Karua, some lava with a Dacite composition of Pleistocene age which is a lava unit of Mt Karua-2, granite rock of Pliocene age which is a Rattebombong intrusion unit, and some parts of Miocene sedimentary rock. According to (Saibi et al., 2008), the high anomaly value can also be interpreted as an indication of a structure that is higher than its surroundings below the surface. This can be indicated by the presence of SE-NW trending structures located in the SE and NE of the study area. While the medium to low anomalies found in the middle to the Western part of the research area are generally located in Miocene andesite rocks which are the lava units of Mt Panusuk, some are in Miocene basalt rocks which are the lava units of Mt Ruppu, some are in basalt, dacite, and andesite of Pleistocene age which is a unit of lava from Mt Karua 1, 2, and 3. This low anomaly can also be indicated by the alteration process in these rocks and is reinforced by the presence of manifestations in the form of fumaroles and hot springs in areas that have low anomalies. 4.2. Horizontal Gradient On the HG map (Fig 4a), the range of HG values ranges from 0 to 0.0085 mGal/m. The anomalous contrast patterns on the HG map are spread throughout the study area. When referring to the structural geological data, the faults located in SW, NE, SE, NW, and the center of the study area shows an anomalous contrast from the HG value with the fault plane at the maximum value of the anomalous contrast. Slicing is carried out on the HG map at the location where the fault is thought to be to ensure and display more accurate results. As previously discussed, there are nine slicing paths with the SVD map at all upward continuation levels. Table 3. Fault characteristics based on derivative analysis. Fault Name Dip Fault Slope Direction Fault Type P01A 64.93º SE Reverse Fault P01B 72.64º NW Normal Fault P01C 66.70º SE Normal Fault P01D 75.38º NW Normal Fault P02A 78.07º NW Normal Fault P02B 56.77º SE Reverse Fault P02C 75.80º NW Normal Fault P02D 64.00º SE Normal Fault P03A 72.88º NW Normal Fault P03B 77.24º SE Reverse Fault P03C 75.61º NW Normal Fault P04A 60.32º SW Normal Fault P04B 74.54º NE Normal Fault P04C 70.93º NE Normal Fault P04D 81.88º SW Reverse Fault P04E 78.26º NE Normal Fault P05A 79.23º NE Reverse Fault P05B 82.50º SW Normal Fault P05C 82.13º NE Normal Fault P06A 76.00º SW Reverse Fault P06B 70.76º NE Reverse Fault P06C 77.51º NW Normal Fault P07A 68.38º SE Reverse Fault P07B 85.84º SE Reverse Fault P07C 60.19º SW Normal Fault P08A 71.05º NW Normal Fault P09A 77.55º NW Normal Fault P09B 84.58º NW Normal Fault 4.3. Multi-Level Second Vertical Derivative The second vertical derivative (SVD) process is carried out on the contours of the gravity anomaly and all contours of the upward continuation of the residual anomaly from a height of 200 m to a height of 2000 m so that 11 SVD maps were sliced together with the HG map. From the 0 m, 200 m, and 2000 m SVD maps (Fig 5), high SVD values are indicated by red color, and low SVD values are indicated by light purple colors, while SVD = 0 is indicated by yellow color (Fig 4b). In the SVD analysis, the value of SVD = 0 indicates the presence of a fault plane and is confirmed by the maximum HG value. Besides that, it is also necessary to pay attention to the maximum and minimum SVD values of the anomalous 190 Putra, U.G. et al./ JGEET Vol 6 No 4/2021 contrast to determine the type of fault. Nine slicing lines were also performed on the anomalous contrasts estimated as subsurface faults (Fig 6). The ML-SVD map is shown in Fig 5. The red color indicates the shallowest structure, while the blue color represents the deepest. On an upward continuous SVD with a height of 200 m, the resulting SVD anomaly contour is almost similar to the SVD anomaly contour at 0 m, where there is still anomalous contrast indicated as a local fault. Meanwhile, at an altitude of 2000 m, the resulting SVD anomaly contour is more regional than the 0 m SVD anomaly contour. This can be seen in the smoother contours and there are only a few anomalous contrast patterns indicated as continuous faults to a depth of 2000 m (middle and north-eastern part of the study area), unlike those found in the 0 m SVD anomaly contour and the SVD contour. Results of upward continuation at an altitude less than 2000 m. From the SVD map at each depth from the upward continuation, the movement of the fault from this shift value SVD = 0 of the SVD map per depth can be seen. From the movement of the value of SVD = 0, the approximate direction and the dip of a fault can be determined. 4.4. Fault Interpretation Based on derivative analysis on nine paths in the Bittuang geothermal area, the correlation between the HG curve and the ML-SVD curve adapted to the structural geological map of the Indonesian Geological Agency (2009) results in an indication of the subsurface structure in the form of 27 faults (Fig 6). Of the 27 faults, nine faults are normal faults with a dip towards NW, two normal faults with a dip towards SE, three normal faults with a dip towards SW, four normal faults with a dip towards NE. As for the rising fault, five faults with a dip to the SE, three faults with a dip to the NE, and one fault with a dip to the SW. Normal faults dominate faults in the Bittuang geothermal area. The result is under the field survey by the Indonesian Geological Agency (2009), which states that the faults in the Bittuang geothermal area are generally normal faults with NE-SW and SE-NW fault directions. From the results of this study, there is also a strike – slip fault with a NE-SW direction that cuts rocks and pre-existing structures. These results are also confirmed by the value of the gravitational anomaly, which does not have anomalous contrasts in the horizontal structure, which is probably along the Eastern to the North-western part of the central part of the Bittuang geothermal area. The gravity method is not suitable for detecting horizontal faults because the gravity method detects the subsurface density contrast on the same horizontal as indicated by the resulting gravitational anomaly contrast value. While the horizontal fault does not have a density contrast in the same horizontal, so the resulting gravity value does not have anomalous contrast. Based on the results obtained from processing satellite gravity data in the Bittuang geothermal area, overall, it conforms with the geological map of the structure of geological research (Indonesian Geological Agency, 2009). Based on the geochemical research, the geothermal manifestations of Bittuang are in two different locations. The first location is in the Balla manifestation group, while the second is in the Cepeng manifestation group. Control faults around the location cause the emergence of geothermal manifestations to the surface. For the Balla manifestation group, the controlling structure of the manifestation is probably the P07A fault located on line 7. The fault is a normal fault with a dip of 68.38º to the NW direction, which continues to a depth of 1400 m from the surface of the residual anomaly. Thus, if later drilling is carried out at a location near the upflow zone, drilling can be carried out in the NW part of the presence of the zone controlling structure (Fault P07A). Meanwhile, information on dip size and fault depth can consider how far and how deep the drilling point is from the control fault location. As for the Cepeng manifestation group, the controlling structure of the manifestation is probably the P09b fault which is on line 9. The fault is a normal fault with a dip of 84.58º to the NW direction, which continues to a depth of 2000 m from the surface of the residual anomaly. Thus, if drilling is carried out in a location near the outflow zone, drilling can be carried out in the NW part of the presence of the zone controlling structure (P09B fault). 5. Conclusion Finally, the ML-SVD method can identify and characterize existing faults in the Bittuang Geothermal prospect area. Based on the analysis of 9 lines, there are 27 faults detected from the gravity data of the GGMplus satellite. Normal faults dominate faults in the Bittuang Geothermal area. The average dip of all faults in the study area is 74.77º with NW, NE, SE, and SW dip directions. The gravity method cannot identify a transform fault that is characterized by the absence of contrasting gravitational anomalies. The P07A fault located on line 7 is indicated as the controlling fault in the Balla manifestation group. The fault is a normal fault with a dip of 68.38º to the NW, which continues to a depth of 1400 m from the surface of the gravity anomaly. Meanwhile, the P09b fault on line 9 is indicated as the controlling fault in the Cepeng manifestation group. The fault is a normal fault with a dip of 84.58º to the northwest, which continues to a depth of 2000 m from the surface of the residual anomaly. Acknowledgments The author would like to thank the Indonesian Geological Agency and enerGIS Indonesia for providing the geological and geochemical data and the Department of Geosciences Universitas Indonesia for research funding. References Abdel, Z.M., 2011. Geothermal Exploration and Numerical Modeling at Gulf of Suez, Egypt. Kyushu University. Blakely, R.J., 1995. Potential Theory in Gravity and Magnetic Applications. Cambridge University Press, New York. Giggenbach, W.F., 1991. Chemical techniques in geothermal exploration, application of geochemistry in geothermal reservoir development. UNITAR/UNDP Cent. Small Energy Resour. Rome 11, 119–133. Hirt, C., Claessens, S., Fecher, T., Kuhn, M., Pail, R., Rexer, M., 2013. New ultrahigh-resolution picture of Earth’s gravity field. Geophys. Res. Lett. 40, 4279–4283. https://doi.org/10.1002/grl.50838 Indonesian Geological Agency, 2009. Integrated Investigations of Geology, Geochemical and Geophysics Bittuang Geothermal Region-Tana Toraja Regency South Sulawesi. Geological Agency, MEMR, Putra, U.G. et al./ JGEET Vol 6 No 4/2021 191 Bandung. Jhanesta, W., Supriyanto, 2021a. Derivative Filter of Gravity Data to Detecting the Structural Control Geothermal Manifestation, Mount Endut GWA, Indonesia, Exploration Geophysics. Jhanesta, W., Supriyanto, 2021b. Application of Multi-Level Second Vertical Derivative (ML-SVD) Method to Identify Subsurface Structure in Mount Endut Geothermal Prospect Area, Indonesia, in: The 2nd Digital Indonesia International Geothermal Convention (DIIGC) 2021. Kusnadi, D., Setiawan, D.I., 2009. Penyelidikan Geokimia Daerah Panasbumi Bittuang Kabupaten Tana Toraja, Provinsi Sulawesi Selatan. Geological Agency, MEMR, Bandung. Nishijima, J., Naritomi, K., 2017. Interpretation of gravity data to delineate underground structure in the Beppu geothermal field, central Kyushu, Japan. J. Hydrol. Reg. Stud. 11, 84–95. https://doi.org/10.1016/j.ejrh.2015.11.022 Rosid, M.S., Naufal, M.A., 2019. Characterization of Geological Structure Based on Multi Level Second Vertical Derivative (ML-SVD) Gravity Data, Journal of Geophysical Research. Saibi, H., Nishijima, J., Hirano, T., Fujimitsu, Y., Ehara, S., 2008. Relation between structure and low- temperature geothermal systems in Fukuoka city, southwestern Japan. Earth, Planets Sp. 60, 821–826. https://doi.org/10.1186/BF03352833 Schon, J.H., 2015. Physical Properties of Rocks - Fundamentals and Principles of Petrophysics. Elsevier. Supriyanto, Noor, T., Suhanto, E., 2017. Analysis of gravity data beneath Endut geothermal prospect using horizontal gradient and Euler deconvolution. AIP Conf. Proc. 1862. https://doi.org/10.1063/1.4991298 Uwiduhaye, J. d. A., Mizunaga, H., Saibi, H., 2018. Geophysical investigation using gravity data in Kinigi geothermal field, northwest Rwanda. J. African Earth Sci. 139, 184–192. https://doi.org/10.1016/j.jafrearsci.2017.12.016 © 2021 Journal of Geoscience, Engineering, Environment and Technology. All rights reserved. This is an open access article distributed under the terms of the CC BY-SA License (http://creativecommons.org/licenses/by- sa/4.0/). http://creativecommons.org/licenses/by-sa/4.0/ http://creativecommons.org/licenses/by-sa/4.0/ http://creativecommons.org/licenses/by-sa/4.0/