Quaternary lava tubes distribution in Jeju Island and their potential deformation risks

Jeju Island, located in the southern sea of the Korean Peninsula, is a volcanic island created by the tertiary and quaternary volcanic eruptions. A group of lava tubes formed between quaternary lava flows constitutes one of the most predominant geological contexts owing to its unique and complex network, for which a total of 178 lava tubes is estimated. 15 As a significant portion of lava caves have not been discovered, the threat caused by lava cave collapse has become one of the major concerns in connection with the recent infrastructure construction in Jeju Island. Considering the risk potential, the overall distribution and collapsing risk of the Jeju lava tube network were investigated in this study. Through spatial analysis, we firstly found that the lava tubes distribution is not correlated with specific geological units. Secondly, the risk of collapse is high especially when there are ongoing artificial constructions around the undisclosed lava tube network. We therefore 20 introduced Interferometric Synthetic Aperture Radar (InSAR) techniques to measure the deformation of the ground surface where lava tube networks distributed underground. InSAR results and the proposed machine learning applications identified that the populations of ground deformations was up to 1-2 mm/year and was inferred to be caused by the instability of the shallow lava cavity. Given that underground cavities could pose serious risks, a detailed physical exploration and threat assessment of potential cave groups is required before intensive anthropogenic construction is developed. 25 https://doi.org/10.5194/nhess-2020-321 Preprint. Discussion started: 17 October 2020 c © Author(s) 2020. CC BY 4.0 License.


Data sets
For this study, it is necessary to carry out the detailed topographic data analyses for tracing of surface deformation as well as 125 subsurface structures. The technology of space geodesy, such as InSAR data interpretation, are being widely used to replace https://doi.org/10.5194/nhess-2020-321 Preprint. Discussion started: 17 October 2020 c Author(s) 2020. CC BY 4.0 License. processing, the optimal data set is Sentinel-1 SAR imagery (Geudtner et al., 2014) as the following bases. First, it is freely available on the public domain, Sentinel data hub. Second, images have been acquired by two SAR-satellite constellations, Sentinel 1A and 1B operating from 2015 and 2016 respectively. Therefore the Sentinel-1 constellation provided the shortest revisiting time (>6 days) among all available InSAR assets. Third, their unique Interferometric Wide-swath mode (IW) operation with Terrain Observation and Progressive Scans SAR (TOPSAR) imaging makes available precise monitoring of 135 the target area with a moderate spatial resolution (20 m in azimuth and a 4 m in range) and 250 km wide swath coverage (Geudtner et al., 2014). Especially its C-band wavelength guarantees relatively minor ionospheric errors carrying longwavelength artifacts compared to L-band InSAR. Since the estimated surface deformation over lava tubes is insignificant compared to the external error components, we were obliged to use InSAR time series analysis to obtain the displacement velocity. In terms of image availability, the image connection geometry in Sentinel-1 InSAR descending modes can be readily 140 constructed covering sufficient period as shown in Fig. 2. It is worthwhile noting that, in order to achieve a comprehensive observation, we employed a strategy combining two different InSAR time-series analysis techniques, one is Persistent Scatterers (PS) to monitor the temporal migration on specific scatterers, such as crossing points of roads and lava tubes.
Another method adopted is Small Baseline Subset (SBAS), which is for the extraction of regional deformation patterns. To achieve sufficient InSAR pair stacks, PS analysis was constructed to cover a two year period with 75 images (Fig. 2(a)), while 145 SBAS network was built on 13 images covering a half year period as shown in Fig. 2(b). The required number of InSAR pairs is usually higher in PS processing. SBAS was mainly introduced to observe regional characteristics compared to PS analysis.
Thus a selection of 13 SBAS images for a half year fits the purpose. Basic characteristics of Sentinel-1 imagery are summarized in Table 1. Another highlight of this study is that a few subsurface morphologies in lava tubes were established by 3D laser scanning to demonstrate their stability. Laser scanning was conducted using FARO Focus3D X330 laser scanner and GPS receiver. Point density of laser scanning was set as 15 mm. A number of targets were pre-installed to be used as a joint point for data merging and also as an inspection point. The targets were installed without interfering with the point cloud acquisition inside the cave.
Once the scanned was finished, the scanned point cloud collected in multiple stations were merged to form one complete model 160 in the FARO SCENE 5.4 data post-processing software. Given the parameters such as the number of points, the number of repetitions and the search range, the registration error value was calculated within ± 2 mm. After then, noise removal and georeferencing were performed based on ground control points to construct a fully co-registered 3D cave model for three lava caves.

Spatial analysis of lava tube distributions
In spite of long historical research works, the distribution pattern of lava tubes still remains unclear. One clue to reveal the pattern is an involvement with low viscosity lava flow so called pahoehoe lava (Ahn and Hwang, 2009). However, as shown in Fig. 1, the spatial extent of pahoehoe lava flow and distribution of lava tubes are not highly correlated. Therefore, we tried 170 to constrain the potential distribution of lava tubes using other evidences. We created the density map of lava tubes (178 in https://doi.org/10.5194/nhess-2020-321 Preprint. Discussion started: 17 October 2020 c Author(s) 2020. CC BY 4.0 License. total) with the weighted values of their lengths and the co-kriging interpolation as shown in Fig. 3(b). The other useful information to spot the subsurface lava cavity is the existence of a geological feature so-called Sumgol, which assigns the collapsing place of lava cavity or the vertical joints or cracks over lava flows (Hamm et al., 2005). It is presumably related to the distribution of undiscovered lava tubes. We digitized the Sumgols and applied the distance transformation with the city 175 block sampling method (Huang and Mitchell, 1994) as shown in Fig. 3(d). The convolution of two spatial distributions was proposed as the representation of subsurface cavities by lava flow (refer Fig. 8(a)). The convolution is corresponding to simple adding of the two materials as we do not apply any weighting. The gridded data sets assigned the distribution of Sumgol and the footprints of lava tubes were employed for the further processing of InSAR data sets as a constraint. The other data sets extracted from the geospatial information were used for the ground truth and validation. The details were described in section 180 4.

InSAR processing
Although InSAR techniques have been effectively used for natural and artificial surface deformation observations, land surface monitoring caused by instability on the subsurface cavity would be a challenge of InSAR techniques due to the expected low https://doi.org/10.5194/nhess-2020-321 Preprint. Discussion started: 17 October 2020 c Author(s) 2020. CC BY 4.0 License.
The PS method exploits strong scatterers in time series interferograms, which are paired by a common master image and sufficient number of slave images (Ferretti et al., 2001;Ferretti et al., 2000). Compared with the PS approach, the SBAS technique (Berardino et al., 2002) determines all the appropriate InSAR pair subsets possessing a small baseline. In this study we employed PS and SBAS to respectively perform a long-term pointwise observation and short-term regional observations. An improved SBAS technique, New-SBAS (NSBAS), was employed to enhance the performance. The algorithmic 195 improvements of time series analysis mainly aim for the densification of reliable scatterers and the capabilities have been proven.
According to a preliminary PS InSAR analysis, it was indicated that the temporal and spatial baseline conditions between SAR images in this study were adequate to address technical challenges with a standard PS algorithm, as shown in Fig. 2(a).
The observations and error component terms to be handled in algorithms of time series analyses is expressed as follows: 200 where ∆φ is the phase difference of interferogram, ∆φ is the phase difference only by target topography, ∆φ is the 205 phase difference by inaccurate orbital information, ∆φ is the phase difference by inaccurate base topography, φ is the other phase noise, ∆φ is the phase difference by atmospheric phase components, and ∆φ is the phase difference caused by delay in the ionosphere.
The core idea of PS algorithm is to discriminate Persistent Scatterers with constant responses for amplitude dispersion and to address the error estimations using iterative non-linear equations. Although PS required a large number of image stacks, the 210 accuracy of InSAR deformation is up to 1 mm/year (Crosetto et al., 2016). On the contrary, PS algorithm is often suffered by the lack of observation density. Thus some PS variants were often introduced. In this study, the target area is full with stable rock scatterers. We therefore are able to employ the conventional PS algorithm to observe the temporal migration with high precision. It was also observed that the density of extracted scatterers was sufficient for further interpretation. Refer to Fig.   4(a) which shows the LOS displacements extracted via the PS time series analysis in the target area. The stable reference point 215 which is as the standard for measuring relative deformations was carefully chosen based on the stability and dispersion of phase coherence and deformation rate. The time series displays of LOS surface deformations demonstrate that the deformation patterns along geological/landcover units were consistent overall. However, there were discrepancies between a few geological units and LOS deformations, particularly along the Manjang, Namhyun and Posunri cave clusters. This implies that there is clear surface migration induced by cave instability. 220 https://doi.org/10.5194/nhess-2020-321 Preprint. Discussion started: 17 October 2020 c Author(s) 2020. CC BY 4.0 License.

225
The outcomes by PS algorithms were further complemented by SBAS technique to manipulate InSAR time series in regional scale. The local deformation was then interpreted together with scatterer behaviours. Since we experienced that SBAS technique is interfered by atmospheric noises, NSBAS technique was employed to improve this performance issue. The https://doi.org/10.5194/nhess-2020-321 Preprint. Discussion started: 17 October 2020 c Author(s) 2020. CC BY 4.0 License. given in Fig. 5. Note that in PS observations, their phase coherences are in high ranges (> 0.75) and the reliability was proved ( Fig. 4(b)). However the phase coherences of NSABS are relatively low except for rocky surfaces and urban areas as shown Fig. 5(b). Hence, we employed NSBAS data as an ancillary data rather than a deformation signal of LTDP. 270 To tackle all problems, we introduced machine learning methods. First, we built the training datasets based on the spatial 290 analyses and the geological/demographic contexts. The established ground truth data sets were then feedforwarded to the training stage of proposed machine learning methods. Subsequently we classified all PS target points using trained machine learning platforms. The extracted points as LTDPs were re-analysed in the comparison of background context and validation data sets.
In Table 2 and Fig. 6, the standards to define training data sets and their locations were shown. As Son (2009) assigned, the 295 road crossing with lava tubes or the cavities on lava flows often caused the failures of the roof of corresponding lava tubes/cavities. We investigated the deformations patterns of PS observations on road-crossing points with Manjang cave and all Sumgols. We choose Manjang cave as the most deformation fragile case because of the shallow structure of the lava cave (see Fig. 9(a)), extensive length and high interaction with tourist activities. On the contrary, Sumgols was originally created by the vertical failures of lava cavities, and the proximity to roads supposedly indicates the instability and deformations. With 300 177 InSAR observations (77 on Manjang cave and 100 on Sumgols) in road crossing areas regions, we discovered that the deformation patterns could be classified as 1) minus LOS migrations (< -1 mm/year) which imply vertical subsidence, so called V-migration patterns; 2) positive LOS migrations (> 0.75 mm/year) which might refer to horizontal deflections, so called Hmigration patterns as consistent uplift deformations are almost impossible especially in rocky surface of lava flow. Thus 50 Vmigrations and 17 H-migrations were assigned as the training data. Note that there might be some minus LOS migrations 305 induced by horizontal creep to LOS direction in SAR sensor but this portion in V-migration should be minor as the major contribution to LOS direction is up-and downward deformation (Hu et al., 2014). Interestingly, around Sumgol, the major of deformations (> 70%) are V-migration, while H-migration in Manjang cave occupied 50% of deformation points. It fits the https://doi.org/10.5194/nhess-2020-321 Preprint. Discussion started: 17 October 2020 c Author(s) 2020. CC BY 4.0 License. hypothesis in which vertical failures in Sumgols mainly induce V-migration patterns, but the large cavity such as Manjang produced H/V-migration pattern together. Therefore, we constructed the training vectors of potential lava tubes combining 310 two V-migrations patterns together. Their behaviours included highly variable signals perhaps involved with seasonal effects as shown in Fig. 7(a) and (b), where changes possibly induced by thermal expansion are noticed in summer time. H-migration patterns were introduced as the secondary indicator of potential instability by lava cavity. The stable PS points with no significant variation were excluded for the input data of machine learning applications. However, we observed that two kinds of deformation patterns might be mis-recognized as the genuine LTDPs. The first one is the deformation induced by regional 315 subsidence, while the other is the structural deformation such as newly built buildings. The regional subsidence is clearly observed in Seowipo sediment perhaps caused by regional condensations (see Fig. 6 and Fig. 1). Therefore the strong minus deformations (< -2 mm/year) in Seowipo sediment were chosen as the training vectors of regional subsidence. The stable behavior is shown in Fig. 7(c) in which the variation is quite different from the Smugol-road-crossing or tube-road-crossing regions. The instability of individual structures was chosen on the opposite side of Jeju city areas where large buildings are 320 populated. Since there are a large number of deformations with random changing patterns, the average behaviour is uniform and is shown in Fig. 7(d). Since there are not enough training points, the pattern of H-movement is highly variable (Fig. 7(e)).

380
Regarding the instability points around Manjang (E) cave, it is corresponding to the exit of Manjang cave which has shallow double cave structure as shown in Fig. 9(a). Groups F1-F3 demonstrates the high density of potential LTDPs and alignment with lava flow directions from cinder cones to coastal line. Groups G1 and G2 are ambiguous as their alignments do not fit with lava flow directions in that area. Although regional subsidence and/or alignment with corresponding lava flow was not observed, a constant attention of the potential LTDP distribution in north-eastern side is still required. Along the north western 385 coastal line, there are some lava tubes which are in collapsing risks due to their shallow cavity and proximity to road-crossing (see the Geamcheon and Gamnamdap caves shown in the laser scanned data respectively given in Fig. 9(b) and (c)). The potential LTDP distributions in the north western coastal line are only limited along the coastal line and irregularly distributed. This could be explained by the deeply-incised caves in inland side to limit the InSAR observation capability in comparison to Geamcheon and Gamnapdap caves. Note that the deep incised caves are located in the inland area from Geamcheon cave. The 390 https://doi.org/10.5194/nhess-2020-321 Preprint. Discussion started: 17 October 2020 c Author(s) 2020. CC BY 4.0 License. issue now being raised is whether these 3D cave structures make deformations as observed in InSAR analyses. Since numerical modelling, for instance the FEM using 3D cave data and estimated tensile stress, is beyond the scope of this study, clues can be found in precedent studies regarding tunnel stability. In particular, Yang and Long (2015) established the relationship between the dimensions of collapsing parts which are defined with L1, L2 and radius R on circularly approximated cavity (see Fig. 9 (d)), stresses and material properties. Thus the displacement induced by collapse will be dominated by the vertical way 395 but mixed with the horizontal deformations which are applied symmetrically at the centre of the approximated circular as shown in Abdellah et al. (2018). Therefore H-migrations detected by PS InSAR mainly contained part of horizontal deformations directed toward the SAR sensor, whereas the detected V-migrations represented mainly vertical deformations combining horizontal deformations away from SAR sensor position. The magnitude of deformations around circular tunnels are very much different according to the applied tensile stress, materials properties (hard/soft rocks and soil) and the dimension 400 of tunnel. If the tunnel depths are limited to very shallow (< 10-20 m), the deformation values calculated in the case studies vary from sub-millimetre depending on the scenarios (Paternesi et al., 2017;Abdellah et al., 2018;Zhang and Li, 2009).
However, the fractures in the weathering wall of the lave tubes and the relatively weak brittle strength of basaltic rocks combined with anthropogenic stresses certainly cause sufficiently high deformations which can be detected by InSAR observations. 405 InSAR LOS migrations can be expressed using the quantities in Fig. 9 (d) and the relationship by Fialko et al. (2001) as (4) where displos is LOS displacement, disph, and dispv are horizontal and vertical displacements, θ represents the incidence angles 410 and α is the heading angle.
Thus, the observed displos can be estimated using the angles in Table 1 and eq. (4) in the case of Fig. 9 (d).
https://doi.org/10.5194/nhess-2020-321 Preprint. Discussion started: 17 October 2020 c Author(s) 2020. CC BY 4.0 License. After all, we identified a group of potential LTDPs using InSAR observation and spatial analyses. To prove the reliability of this approach, we employed road-crossing survey by Son (2016) to measure the known collapses of lava tubes. Among 27 ongoing collapses of lava tubes in this survey, 14 caves have LTDP within the buffer zones built by their length and 7 of them have both H-V migrations (referred to Fig. 10). 430 https://doi.org/10.5194/nhess-2020-321 Preprint. Discussion started: 17 October 2020 c Author(s) 2020. CC BY 4.0 License.