Recent progress of volcano deformation studies

A lecture on volcano deformation studies given in the Asia Volcano Consortium Field Camp at Academia Sinica.

Transcript

1. Recent progress of volcano deformation studies Yosuke Aoki Earthquake Research Institute, The University of Tokyo Email: yaoki@eri.u-tokyo.ac.jp 29 October 2019 ACV-FC4 Academia Sinica, Taipei, Taiwan
2. Schedule this afternoon 1:00-2:00pm (could be longer) “Recent progress of volcano deformation studies” 2:00-2:30pm Break 2:30-4:00pm (could be shorter) Exercise: “Theory and application of Synthetic Aperture Radar” 4:00-4:30pm “Geodetic analysis of Tatun Volcano” by Masayuki Murase 4:30-5:10pm Poster
3. Why volcano deformation matters? Figure 1. (left) Conceptual diagram (not to scale) of magmatic–hydrothermal interactions (e.g., fluid FOURNIER AND CHARDOT: HYDROTHERMAL DEFORMATION INVERSION B1 8 Fournier & Chardot (JGR, 2012) ü Where and what shape is (are) pressure source(s)? ü How do they move?
4. Development of Volcano Geodesy Development of observational techniques (how to measure) - Conventional methods (leveling, tiltmeters, strainmeters) are still important. - Space geodetic methods (GNSS, SAR) dominate the current volcano geodesy. - Modern techniques (Ground-based SAR, Airborne SAR, LiDAR, SfM) might change the world of volcano geodesy in the future. Development of modeling methodology (how to model) - Simple sources (sphere, dike) embedded in a homogeneous, elastic, and isotropic halfspace. - Numerical modeling (e.g., FEM) and analog modeling incorporating material heterogeneity, topography, or complex source geometry. - How to invert for the source properties with a complex problem setting?
5. Volcano geodesy with conventional techniques Iguchi (Bull. Volcanol. Soc. Jpn., 2013) Mogi (Bull. Earthq. Res. Inst., 1958) ü Less temporal resolution than modern data ü Still important because of longer-term measurements
6. Emergence of space geodetic techniques Aoki et al. (Science, 1999) ü More temporal resolution and temporal stability with GPS. ü Tilt data make an important role in constraining the magma transport during a seismic swarm.
7. Where is the magma reservoir? Eruptions in Aug. 2008 and Feb. 2009 Enhanced seismicity Summer 2008- Shallow inflation Summer-winter 2008 Nagaoka, Nishida, Aoki et al. (EPSL, 2012) Aoki et al. (Geol. Soc. Lond. Spec. Publ., 2013) Asama Volcano
8. Emergence of Synthetic Aperture Radar Amelung et al. (Nature, 2000) ü Spatial resolution unattained by ground-based monitoring with SAR. ü No need to install equipments on the ground. ü Spatial resolution only down to few to few tens of days depending on the recurrence time of the satellite.
9. Volcano deformation of various origins A Concentric C Flow Subsidence D Landslides G Lateral Feeding System F Vertically Stacked E Cooling Intrusion B Dyke Intrusion Emergence of InSAR revealed that volcano deformation is more complicated than we thought. Biggs & Pritchard (Elements, 2017)
10. Heuristic approach of volcano deformation Pritchard and Simons (Nature, 2002) Chaussard, Amelung, and Aoki (JGR Solid Earth, 2013)
11. Global compilation of deformation sources ween the centre of deformation and the nearest edifice catalogued in the VOTW database are shown with respect to inferred ü 24 % of volcano deformation takes place >5 km from the volcanic center. ü Majority of the deformation source is <10 km. Biggs & Pritchard (Elements, 2017)
12. Global compilation of deformation rate and duration A compilation of ground- based and space geodetic techniques show that slower deformation lasts longer. lava/pyroclastic flow, or pressurizing ‘magma re all of which have distinct patterns of ground d ment viewed with InSAR (FIG. 1A-E). Complexit depending on the direction the satellite is looking to the ground displacement (e.g. Dzurisin et al. 20 when multiple processes occur at nearly the sa (FIGS. 1F, 1G). In this section, we describe some com ties between deforming volcanoes and global pat the types of processes that produce ground defor The Classic Volcano Deformation Cycle Prior to an eruption, according to the classic m the ‘volcano deformation cycle’, magma gradually a magma chamber directly beneath the volcani until a threshold is reached at which point the c ruptures and an eruption rapidly empties and defl chamber (e.g. Dzurisin et al. 2006). The inflatio causes uplift of the ground surface and large num small earthquakes (volcano-tectonic seismicity eruption is accompanied by rapid subsidence (FIG. first approximation, this pattern of co-eruptive su 10−1 100 101 102 103 104 105 10−1 100 101 102 103 104 105 106 Duration (days) Magnitude (mm/y) Eruption No eruption 1 year 10 years 100 years Biggs & Pritchard (Elements, 2017)
13. Global compilation of deformation and eruption Biggs et al. (Nat. Comm., 2014) Red: deformed Black: not deformed Blue: eruption ing field measurements, the variation in such diagnostic test istics associated with globally mapped volcanological and onic factors could provide a solid footing for hazard analysis. quantify the causal or temporal relationships between the extent, rates, timing, magnitude and duration of deformation on shorter timescales. Global studies of volcano deformation have potential to be incorporated into strategic hazard assessments25, particularly in regions with short historical records15. Methods Volcano catalogue. There are 1390 named subaerial Holocene volcanoes13 of which 161 have InSAR-based reports of deformation to date and 620 have reported InSAR observations (Fig. 1). Table 1 and Supplementary Table 1 list systematic studies that include discussion of null results: together they cover 4500 volcanoes in the Andes, Central America, Alaska, Africa, Indonesia, Iceland and the Galapagos. The volcanoes in Iceland are not covered by a single systematic study but 85% of them have been included in separate studies of volcanic, seismic, cryospheric or geothermal processes26–29 and those in Galapagos are sufficiently close together to be covered by a single satellite frame. InSAR studies outside these regions tend to focus on individual events, such that the proportion of volcanoes that erupted (32%) or deformed (59%) is significantly higher than the corresponding values of 12% and 17% for the systematic studies. Ground-based geodetic networks provide valuable information at higher temporal resolution and over a longer time period than is accessible using InSAR. However, due to the high logistical and financial overheads, they exist only for a small number of volcanoes, often those that are known to be deforming or erupting. Due to the bias away from ‘null results’ associated with individual studies and ground-based networks, we base our subsequent analysis on systematic InSAR studies alone (Table 1). Due to the high spatial resolution and coverage, InSAR has been very successful in detecting deformation within volcanically active regions. Sometimes, however, this deformation is up to 25 km from any catalogued volcano summit13. In these cases, the deformation has been attributed to the closest listed volcano. Supplementary Table 6 lists deformation attributed to a volcano not in the GVP list 25 29 9 135 Non-deformed Deformed Non-Erupted Erupted Systematic Coverage True positive False positive False negative True negative DE DE DE DE re 2 | Contingency table linking volcanoes that deformed and pted. The table reports the number of occurrences in each category of the 198 systematically observed volcanoes over 18 years (Table 1, plementary Table 1). A volcano that ‘deformed’, D, is one where at least period of deformation has been observed with InSAR, while ‘not rmed’,  D, means InSAR measurements were made, but that no rmation was reported. ‘Erupted’, E, and ‘not erupted’,  E, volcanoes are e that erupted or not13. See Supplementary Tables 2–5 for details of vidual volcanoes. O 0 1 Proportion 3 Figure 5 | Effect of obse window length is shown volcanoes and the value the 198 volcanoes for w the PPV to increase to on over which decisions are Eruption +5 +10 +15 –5 –10 –15 Years before Years after Refs Alu, 2008 Dalaffilla, 2008 1997 2008 1995 Fernandina, 2005 2009 Augustine, 2006 Hekla, 2000 2003 Llaima, 2007 2008 Sierra Negra, 2005 1998 2008 Alcedo, 1993 Akutan, 1992 Westdahl, 1991 Kanaga, 1994 Seguam, 1993 Gareloi, 1996 Fourpeaked, 2006 Fogo, 1995 Eyjafjallajokull, 2010 El Hierro, 2012 Dabbahu, 2005 Manda Harraro, 2007- Amukta, 1996 2006 1998 1992 1994 1995 2000 2000 2002 Tinguiririca, 1994 Chillan, 2003 Chaiten, 2008 Kizimen, 2010 Okmok, Cerro Azul, Atka, Copahue, Sabancaya, 5 45 46 47 48 49 21 51 52 53 54-56 60 57 59 60 61 62 64 63 40 40 40 40 40 40 34 Maule earthquake 31,34 NATURE COMMUNICATIONS | DOI: 10.1038/ncomms4471
14. Temporal evolution? Tiltmeters installed on the crater rim at volcano (Montserrat) in 1996–1997 det tion patterns correlated with seismicity (Voight et al. 1999). The debate continue these observations should be interpreted a of shallow reservoirs, shear along the cond passage of volatiles (e.g. Nishimura 2009). Lava lakes are essentially magma-filled co the surface, and, in a simple conceptual mo the reservoir pressure would result in chan rather than surface deformation. Many of lakes happen to be in remote or inaccessib the best studied is undoubtedly Halemau where Patrick et al. (2015) showed that d lake levels occur simultaneously – eviden conduit linking the lava lake to the magm Restless Calderas Caldera systems have long repose period large eruptions, but do not remain qui FIGURE 3 Temporal patterns of volcano deformation (schematic); stars represent eruptions. (A) Classic A D B E C F Biggs & Pritchard (Elements, 2017) Temporal evolution of volcano deformation is variable.
15. A tale of two failed eruptions Hotta et al. (Earth Planet. Space, 2016) During the period with the highest deformation rate from 10:27 to 11:54, the seismicity drastically increased, including the largest earthquake at 10:47 (M2.3). After 11:54, the seismicity decreased except for the felt earth- quakes at 14:39, 14:46 and 14:52, after which, only 73 earthquakes occurred on August 16. Such large deformation, as indicated by the change in the tilt and strain, was also detected by GNSS. Temporal changes in daily positions at some sites relative to SVOG, in addition to changes in the tilt and strain at the HVOT, during the period from July 1 to September 1, 2015, are shown in Fig. 5. e positions of the GNSS sites suddenly shifted on August 15, 2015. Before August 14, no signifi- cant deformation was detected. During the period from August 14 to 16, more than half of the stations moved significantly: KURG, located in eastern Sakurajima, moved 6 cm eastward; ARIG, located in southeastern Sakurajima, moved 3 cm toward the south; and HARG 2.7 km northwest of the summit crater of Minami-dake moved 2 and 4 cm toward the north and west, respec- tively. Conversely, north to northeastern (FUTG, KMNG) and southwestern (SBTG) parts of Sakurajima showed Fig. 3 Seismicity and ground deformation during the period from Fig. 4 Hourly number of VT earthquakes during the period from 5:00 to 23:59 on August 15, 2015 Aoki et al. (Science, 1999) 2015 Sakurajima 1997 Izu Peninsula
16. Temporal evolution of volcano deformation from InSAR Aoki and Sidiq (JVGR, 2014) Large number of image acquisitions enables us to retrieve time-varying deformation from SAR images. generated from JERS images contain large topographic residuals near the vents of 2000 eruption. The areas with large topographic residuals might cause wrong phase unwrapping so are therefore masked (the white rectangle in Figure 5b). We observe two localized LOS subsidence from the ALOS‐1 images: one is a more widespread region around the NC crater, and the other one is a more localized deformation around the KC crater (Figure 1). Both the ascending and descending images detected significant deformation around the NC crater, while only the descending images detected significant deformation around the KC crater. Figure 5. (a) The digital elevation model used in the interferometric synthetic aperture radar processing. The black dashed polygons represent vents of the previous eruptions. P1P1′, P2P2′, and P3P3′ are profiles across the 2000, 1977, and 1943 eruption sites, respectively. (b) Mean line‐of‐sight velocities derived from JERS images using the stacking method. The white rectangle indicates the masked region due to large topographic errors. (c, d) Mean LOS velocities derived from the small baseline subset processing of ascending and descending ALOS‐1 images. (e, f) Mean line‐of‐sight velocities obtained from the SBAS processing of ascending and descending ALOS‐2 images. The time span in each sub- figure indicates the observation period. The yellow triangles represent the dome summits that are the same as those in Figure 1. NC = Nishiyama Crater; KC = Konpirayama Crater; KU = Ko‐Usu; OU = O‐Usu; US = Usu‐Shinzan; SS = Showa‐Shinzan. 10.1029/2018JB016729 al of Geophysical Research: Solid Earth 7.2. Comparison of Source Volumes With the Other Geophysical Studies The volumes of heat sources inferred from our thermoelastic modeling for the 1977 and 1943 eruption sites are comparable to the previous estimates. There is, however, a significant discrepancy for the 2000 eruption site in that our estimate is much smaller than that inferred from the coeruptive model. Here we explain how to reconcile this discrepancy. Seismic studies have revealed that a magma chamber at a depth of about 5,000 m below the summit crater is responsible for the activity (Yamamoto et al., 2002). They also revealed that the magma from this level migrated upward to the northwest toward the 2000 eruption site. Calculations show that the ground vertical subsidence in response to thermoelastic deflation due to a heat source of 100 × 106 m3 at a depth of 5,000 m becomes less than 1 mm/year after 6 years of emplacement. Therefore, the inferred volume of the heat source associated with the 2000 eruption solely corresponds to the magma transported into the shallow depths, while the value given by Miura and Niida (2002) corresponds to the total volume change involved during the 2000 eruption. Terada et al. (2008) also pointed out this inconsistency, and they estimated a heat source volume of about 8 × 106 m3 that is quite similar to our estimate from the observations of heat dis- charge rate near the 2000 eruption vents. Although our estimates of intruded volumes for the 1977 and 1943 eruption sites are comparable to the previous estimates, we cannot rule out that there were no deep intrusions during these two unrests because of the limited geodetic measurements at the time. The discrepancies in parameters of our estimates from other geophysical studies may partly result from the limitations of our thermoelastic model. First, mechanisms other than thermoelastic deflation may also con- tribute to the observed subsidence, although they are not the controlling factor as discussed in section 5. Second, our modeling builds on several approximations to simplify the problem. We assume spherical heat sources embedded in a homogeneous, isotopic, and elastic medium. Although a medium with depth‐varying or three‐dimensional variations of material properties and magma bodies with a more complex shape must be a better approximation, we rely on the simple approximation to keep the problem simpler. Also, the ther- moelastic modeling employed in this study (Furuya, 2005) does not involve configuration of temperature of rocks surrounding the magma body, which may also influence the modeling results. Figure 15. Same as Figure 13 but with the 1943 eruption site. The gray triangles and filled cycles denote the leveling and Interferometric Synthetic Aperture Radar measurements, respectively. 10.1029/2018JB016729 Journal of Geophysical Research: Solid Earth Wang and Aoki (JGR Solid Earth, 2019)
17. Classic methods of volcano deformation modeling ü Analytical solutions are available for many kinds of geometry of pressure sources embedded in an elastic, homogeneous, and isotropic halfspace. ü Availability of analytical solutions is a big advantage to be included in inverse problems. Bonaccorso and Davis (JGR, 1999) Segall (2010) Nishimura (JVGR, 2009) Yang et al. (JGR, 1988) Fialko et al. (GJI, 2001)
18. Adding complexity Medium - Topography - Heterogeneous elastic properties - Visco(plasto)elasticity Source - Arbitrary geometry Methods - Analytic - Numerical simulation (e.g., FEM) - Lab experiments
19. Deformation modeling incorporating topography Classic methods for correcting topography ü Williams and Wadge (JGR, 2000) gave an analytical solution of deformation field incorporating arbitrary topography as long as (H/L)2<<1, where H and L represent the scale of topography change and horizontal scale, respectively. ü Topography needs to be taken into account (only) when the depth of the pressure source is comparable with typical horizontal scale of topography change.
20. Deformation modeling incorporating heterogeneous elastic constants üIf the elastic constants are function of depth only, the surface deformation field can be analytically derived by the Thomson-Haskell method widely applied in seismology with a zero frequency limit (Zhu and Rivera, GJI, 2002). üIf the elastic constants vary in horizontal direction as well, the (approximated) analytic solution still exists as long as perturbation is small (Du et al., GRL, 1997; Cervelli et al., JGR, 1999). üIgnoring heterogeneous elastic constants in vertical direction underestimates the depth and volume change of the pressure source (Manconi et al., GJI, 2007; Long and Grosfils, JVGR, 2009; de Zeeuw-van Dalfsen et al., JVGR, 2012). üIt is important to consider (at least) vertical variation of elastic constants in modeling deformation field.
21. Deformation modeling incorporating viscoelasticity “Effective” size of the spherical reservoir becomes progressively larger over time (e.g., Dragoni and Magnanensi, PEPI, 1989) A viscoelastic halfspace overlained by an elastic layer can explain a deflation after a magma injection which is observed in many active volcanoes.
22. Modeling the observation by thermoelastic deformation of intruded lava dome Wang & Aoki (JGR Solid Earth, 2019) V d Sea level Intruded magma body Surface Thermal diffusion Temperature Time elapse High Low
23. Deformation modeling with numerical methods Finite Element Method Masterlark et al. (JGR, 2012) Boundary Element Method Maccaferri, Rivalta, Passarelli, and Aoki (EPSL, 2016) Numerical methods are capable of incorporating complex source geometry and material properties.
24. Insights from analog modeling Trippanera et al. (JGR Solid Earth, 2015) Lab experiments can capture features that sometimes cannot be found by numerical simulations.
25. Combining physical models with observations Anderson and Segall (JGR Solid Earth, 2013)
26. Observational techniques in the future Ground-based SAR Airborne SAR
27. Ground-based InSAR at Stromboli derived from the integration of three different sensors (Hu et al., 2014; Bardi et al., 2016). However, in the Stromboli case, the integration of the GBInSAR and CSK-SqueeSAR™ displacement vectors allows the estima- tion of the displacement in two dimensions (Fig. 8). Displacement vectors measured by the GBInSAR devices are directed towards the sensor, whereas the CSK-SqueeSAR™ displacement vectors are directed away from the satellites, and the vector reconstructed in two dimensions suggests a movement along the slope (Fig. 8). In the of slow displacement rates, it would account for a general creep o SdF, as suggested by Intrieri et al. (2013) and Nolesini et al. (2013 GBInSAR and CSK-SqueeSAR™ displacement time series of the area are compared in Fig. 8. In the case of the CSK-SqueeSAR™ disp ment time series, areal averaging was applied (Fig. 5b). Despite difference in the LOS direction, the two time series are fully compar Fig. 7. a) Cumulative LOS displacement recorded by the GBInSAR during the interval of January 1st, 2010–December 14th, 2014; b) cumulative LOS displacement recorded by the GB during the period of August 7th, 2014–August 9th, 2014. The cumulative displacement map is strongly affected by temporal decorrelation due to rapid movements related to bo sliding of the NEC debris talus and the fast lava flows; c) GBInSAR interferogram formed using a pair of images recorded between 12:25 GMT and 14:38 GMT on September 10th, (first pulse recorded by the system). The interferogram is affected by phase wrapping. In the picture, the different sectors (crater terrace, debris talus, SdF1, and SdF2) are show well as the position of the time-series analysis point. The images are overlain on a 0.5 × 0.5 m DEM-generated hill-shade image. 106 F. Di Traglia et al. / Geomorphology 300 (2018) 95–112 and lava overflows caused frequent gravel flows on the steep SdF slope (Di Traglia et al., 2014b; Calvari et al., 2016). 3. Methodological approach 3.1. Morphological analysis A detailed morphometric analysis of the SdF was carried out using a high-resolution DEM with a spatial resolution of 50 cm (see Salvatici et al., 2016 for the DEM description). To derive the SdF terrain rough- ness map, any holes in the DEM were initially filled (“no data” pixels), and then, the Topographic Position Index (TPI; Guisan et al., 1999; Gallant and Wilson, 2000; Jenness, 2006) was calculated. TPI measures the relative topographic position of the central point as the difference between its elevation and the mean elevation within a defined neighbourhood (3 pixel radius) (De Reu et al., 2013). The TPI was calculated using the focal operators and raster calculator of ESRI's ArcGIS 10.0 using the following equation: TPI ¼ ELE−minELE ð Þ= maxELE−minELE ð Þ ð1Þ where ELE is the DEM-derived elevation and minELE and maxELE are the minimum and maximum elevation values, respectively. Both minELE and maxELE were calculated on a 3 × 3 pixel kernel window. Positive TPI values imply that the central point is located at a higher elevation than its average surroundings, while negative values denote a height lower than the average. The TPI depends on both elevation differences and the predetermined analysis radius/window (De Reu et al., 2013). 3.2. Detecting lithological and morphological changes by means of spaceborne SAR amplitude images A dataset of 85 COSMO-SkyMed SAR (CSK) images acquired in descending orbit between February 22nd, 2010, and December 18th, 2014, was adopted for this study. The CSK images were collected in STRIPMAP - HIMAGE mode, achieving medium resolution, wide swath imaging and VV single polarization. The swath extension is ≥40 km, with a spatial resolution of 3 × 3 m single look, which is optimal for land-cover change detection (Covello et al., 2010; Bignami et al., 2013). The images were exploited as standard products, processing Level 1A - Single-Look Complex (SLC) slant products, with RAW data Fig. 1. a) Location of the island of Stromboli. b) The island of Stromboli. The Sciara del Fuoco (SdF), crater sectors and the location of the GBInSAR device are shown. 97 F. Di Traglia et al. / Geomorphology 300 (2018) 95–112 Di Traglia et al. (Sci. Rep., 2015; Geomorphology, 2018) .com/scientificreports/
28. Airborne SAR in Kilauea While a SAR satellite looks at the target from east or west because it takes a polar orbit, airborne SAR is capable of looking at the target from any direction. Lundgren et al. (JGR Solid Earth, 2013)
29. Detecting volcanic plumes with GNSS ummary, we propose that the eruption pulses T3 and T4 each ely ejected a dense ash plume that intersected the VGKR-to- th, respectively causing the observed sub-phases Ia and Ib in ata. II avalanche, and the later NW-directed explosion (T3) with the associat- ed pyroclastic density current (PDC; see Lube et al., 2014–in this issue). To discriminate between these two potential sources, one must first as- sess what physical processes induced the increase in absolute residuals. The computed residuals are an expression of a delay in radio wave transmission between the satellites and the GPS antenna. Such delay can be caused by changes in atmospheric conditions along the path of PRN03 PRN13 PRN03 PRN13 35 40 45 50 55 60 65 db−Hz 11:45:00 11:30:00 L2 SNR data for stations VGKR & VGOT 35 40 45 50 db−Hz 11:52:00 11:52:00 11:52:00 11:52:00 11:52:00 PRN03 PRN13 eruption onset T1 T4 T3 Ia Ib VGKR L2 SNR data ata. Top: 40 min worth of SNR data for satellites 3 and 13 at the two GPS stations VGKR (black) and VGOT (blue). Note the drop in SNR data for PRN03 at VGKR only. The box ~3-min time window containing the eruption and is zoomed in the bottom panel. SNR data of PRN03 for VGOT, PRN13 for VGKR and VGOT are offset by 4, 10 and 14 db-Hz to add clarity. Bottom: zoom of the SNR data for PRN03 and PRN13 at VGKR for an ~3-min time window encompassing the eruption. Note the two distinct negative anomalies Ia and Ib) in SNR data for PRN03. The onset of the eruption and of each main pulse is indicated as a vertical marker (see text for more information). 393 N. Fournier, A.D. Jolly / Journal of Volcanology and Geothermal Research 286 (2014) 387–396 5 min. Large particles degrade the strength of GNSS signals. GNSS can infer the size distribution of ejecta? 2012 Tongariro (Fournier & Jolly, JVGR, 2014) 2015 Kuchinoerabu Aoki & Larson (in prep.)
30. Modeling deformation associated with hydrothermal activity responsible for ground deformation (Figure 1)? Ground sur- face deformation caused by hydrothermal activity is driven by of this work, and draw some conclusions about how our results can be used in real–life case studies. Figure 1. (left) Conceptual diagram (not to scale) of magmatic–hydrothermal interactions (e.g., fluid migration) and resulting ground deformation due to hydrothermal unrest. (right) Geophysical observables (e.g., ground deformation) in this case and the information derived from their inversion (e.g., location, size and overpressure of a simplified pressure–source). FOURNIER AND CHARDOT: HYDROTHERMAL DEFORMATION INVERSION B11208 B11208 Difficult to identify pressure source. An elastic modeling is enough? Fournier & Chardot (JGR, 2012)
31. Deformation associated with phreatic eruption: 2014 Ontake post eruptive deformation matrix as suggested in Morishita et al. (2016), whereas non-diagonal components were zero. A coherence value computed over a small spatial window around the con- sidered pixel was used in the variance estimation process. of the deflation (Fig. 3). The deflation began the day after the eruption and continues as of July 2017. The largest dLOS velocity values in coherent areas were 21 cm/year for path 126 and 12 cm/year for path 20, respectively. The temporal profile of the dLOS for the - 2015/7/27 d Path 27: ErupƟon -2017/5/28 a Path 126: 2014/10/3-2017/7/21 e Path 121: 2014/10/6 2014/9/28-2015/7/19 c Path 20: 2014/10/5 b Path 19: 2014/9/30-2016/9/13 f -12 cm +12 cm 0 2km Away Toward Fig. 3 Results of SAR interferometry for the area in Fig. 1 inset. a–e are interferograms depicting the stacked dLOS. The path and period for the stacking are shown at the top of each interferogram. f Time series of the dLOS at fixed pixels as indicated with black circles in (a–c). Positive dLOS values correspond to increases in the slant range between the satellite antenna and the ground target. Arrows represent the flight and viewing direction of the ALOS-2. Black dots are eruptive vents for the 2014 eruption source was parallel to the vertical axis (i.e., a non-dipping ellipsoid) and that the spheroid was axisymmetric around We tested our bes pendent data set. e west east -20 down (cm) up EW UD Observation FE Model b a c d 2 km 2 km 2 km 2 km Fig. 5 Results of FEM modeling for the area in Fig. 1 inset. a, b Observations of the EW and UD compo Fig. 4a, c, respectively. c, d Synthesized distributions for the corresponding components based on the combinations of the free parameters, that is, the aspect ratio of the spheroidal source (abscissa: ratio o shape schematics) and the depth (ordinate). In the left half of the chart, the aspect ratio is represented axis and the shape is characterized as prolate (vertically elongated), whereas in the right half of the ch vertical axis and the shape is oblate. The residuals are colored in grayscale. A red star indicates the bes To improve the fit, we then built a finite element model (FEM). We employed the Salome-Meca software package (https ://www.code-aster .org) to (1) take topo- graphic effects into account and (2) to allow flexibility for the source shape. For the topographic data, we used a 0.4-arc-sec DEM generated by GSI. We set the resolution as 20 m in the 4 km × 4 km region around the Jigokudani Valley and 2 km in the outer margin region. e size of the entire model was 26 km × 26 km × 8 km. We placed a spheroidal source with horizontal position fixed to the best-fit position of the Mogi source. e spatial pattern of the deformation appeared axisymmetric. For simplic- ity, we assumed that one of the axes of the spheroidal source was parallel to the vertical axis (i.e., a non-dipping ellipsoid) and that the spheroid was axisymmetric around tion volume between October 2014 and July 2015 was 3.5 × 105 m3, and the cumulative volume change through 2017 was 7 × 105 m3 assuming the same temporal profile as that of dLOS change in path 126 (Fig. 3f). Our mod- eled source was emplaced at a shallower depth than any deformation and seismic sources reported by other authors (Fig. 6) such as the tilt source at 1100 m depth prior to the 2014 eruption (Maeda et al. 2017), the infla- tion source at 900 m and 3 km depth between 2005 and 2007 (Takagi and Onizawa 2016), the subsidence source at 6 km depth revealed by leveling after the 2014 eruption (Murase et al. 2016) and the VT hypocenters, distributed around 3 km depth (Kato et al. 2015). We tested our best-fit FEM source using another inde- pendent data set. e Japanese Metrological Agency e EW UD Observation E Model b a c d 2 km 2 km Narita & Murakami (EPS, 2018) 0.75 Mm3 of deflation in the first 3 years after the eruption
32. Deformation associated with phreatic eruption: 2015 Hakone Doke et al. (EPS, 2018) in the Owakudani area were detected by ground-based SAR observation almost at the same time (Doke et al. 2015b). Mannen et al. (2015a) reported that volcanic ashfall started at around 0:30 p.m. (JST) on June 29. e observes surface displacements by multi-mode, right- and left-looking, multi-off nadir angles with a 14-day revisit interval. ALOS-2/PALSAR-2 observed the Hakone Volcano under five different orbits during the period of Fig. 3 Temporal change of local swelling at Owakudani (modified from Doke et al. 2015a). Interferograms are obtained from InSAR pairs in the same orbit with Pair I (Path 18). Color shows wrapped phase difference between 0 and +2π radian. Phase difference of 2π radian corresponds to a slant-range change of 11.9 cm. Contour intervals are 25 m in height Table 3 Open crack and sill (deflation source at depth) parameters estimated from the second inversion a Altitudes are height above sea level of the crack top edges, with ellipsoidal heights in parentheses b East and North coordinates are in WGS84 and refer to the upper left corner of cracks c Parameter was fixed Length (m) Width (m) Altitudea (m) Strike (°) Dip (°) Eastb (°) Northb (°) Opening (cm) 1192.1 298.6 827.6 (868.1) 323.83 88.52 139.030707 35.233585 31.86 261.4 304.1 224.9 (265.4) 0.00 0c 139.027414 35.235381 – 62.10 Fig. 7 Result of the inversion estimating the open crack and deflation sill model. a Modeled surface displacements, and b residuals between model and observation for Pair I. c Modeled surface displacements, and d residuals between model and observation for Pair II. B and C indicate displace- ment areas mentioned in the text. KM and SU show the summits of Kamiyama and Sounzan, respectively. Offset parameters were 1.95 mm for Pair I and –9.20 mm for Pair II. The red line and blue rectangle show the estimated open crack and deflation sill, the parameters of which are listed in Table 3 The observation is modeled by a dike intrusion. Page 12 of 14 ke et al. Earth, Planets and Space (2018) 70:63 ig. 10 Conceptual diagram of the 2015 volcanic activity of Hakone Volcano inferred from the data modeling of InSAR images. a The reservoir of
33. Beyond elastic modeling with a generic model. All models are thus deliberately kept simple but bear important characteristics common to most known hydrothermal systems. [7] The hydrothermal fluid–flow modeling strategy is based approximated as a point source, which is consistent with the injection of magmatic fluids at the base of a hydrothermal system through a fracture, as the source of hydrothermal unrest. [11] 4. The fluids are multicomponent (H2 O and CO2 ) and multiphase (liquid and gas). [12] 5. The fluids are injected at a constant rate throughout the simulation period. [13] The base model (1) is a reproduction of one of the models by Hutnak et al. [2009] and depicts hydrothermal system dynamics for a large caldera system (50 Â 5 km, radius and depth respectively). Because volcanoes and associated hydrothermal systems’ dimensions span across several orders of magnitude (i.e., O(102–104) m), we tested the influence of the size of the modeled system on the simulations results. Model 2 (intermediate, 5 Â 1.7 km) was designed with smaller dimensions and same injection characteristics (i.e., injection rate and temperature). Model 3 (intermediate, 5 Â 1.7 km) has the same dimensions as model 2 but a lower injection rate to test the effect of injection rate on our results. Finally, model 4 (small, 5 Â 0.5 km) has smaller dimensions and lower injec- tion temperature. Dimensions and injection characteristics for the three models are summarized in Table 1. 3.1. Heat and Fluid–Flow [14] Heat and fluid–flow modeling was conducted using the TOUGH2 code [Pruess, 2004] with the EOS2 module. It allows simulations of nonisothermal flows of multicompo- nent, multiphase fluids in multidimensional porous and fractured media. The reader is directed to Christenson et al. [2010], Hurwitz et al. [2007], Hutnak et al. [2009], Rinaldi et al. [2010], Todesco [2009], and Todesco et al. [2010] and references therein for extensive details of the TOUGH2 modeling for volcanic hydrothermal systems. [15] Boundary conditions are imposed at the ground surface (fixed temperature and pressure of 10C and 0.1 MPa respec- tively in our case), outer edge and at the center of the model (symmetry axis). The lower boundary represents the base of the hydrothermal system: it is impermeable, apart from the injec- tion point where fluids are injected into the domain, and has constant basal heat flux (0.1 W mÀ2 in our case). Initial con- ditions in the domain are hydrostatic pressure and linear vertical thermal gradient between top and bottom boundaries, as in Hurwitz et al. [2007] and Hutnak et al. [2009]. Figure 2. (top) Flowchart showing the different steps used for simulating hydrothermal models and calculating the deformation that they induce and (bottom) the methods for retrieving each model characteristics from the synthetic sur- face displacements. FOURNIER AND CHARDOT: HYDROTHERMAL DEFORMATION INVERSION B11208 B11208 [2007] and Hutnak et al. [2009], pore–pressure increases sig- nificantly at – or near – the onset of the simulation, while the temperature increase remains limited. Through time, pore– pressure generally decreases near the injection point and gets redistributed as fluids migrate more evenly throughout ground deformation, we computed for each model the ground deformation due to the pore–pressure effect only. This was achieved using an unrealistic low thermal expansion coeffi- cient (a in equation (5)) in BIOT2, hence effectively inhibiting thermal expansion in the models. We compared the results of Figure 3. Example of temporal evolution of (top) pore pressure and (bottom) temperature in the TOUGH2 models (e.g., caldera model at early (1’s years) and later stages (100–1000’s years); note that this plot is a close–up of the first 3 km from the injection point, while the model’s horizontal dimension is 50 km). Through time, pore pressure first accumulates near the injection point and then get redistributed into the model’s domain, with a size reduction of the pressurized zone near the injection point. Conversely, the zone of high temperature expands through time, with the focus remaining the injection point. FOURNIER AND CHARDOT: HYDROTHERMAL DEFORMATION INVERSION B11208 B11208 7.2. Comparison of Source Volumes With the Other Geophysical Studies The volumes of heat sources inferred from our thermoelastic modeling for the 1977 and 1943 eruption sites are comparable to the previous estimates. There is, however, a significant discrepancy for the 2000 eruption site in that our estimate is much smaller than that inferred from the coeruptive model. Here we explain how to reconcile this discrepancy. Seismic studies have revealed that a magma chamber at a depth of about 5,000 m below the summit crater is responsible for the activity (Yamamoto et al., 2002). They also revealed that the magma from this level migrated upward to the northwest toward the 2000 eruption site. Calculations show that the ground vertical subsidence in response to thermoelastic deflation due to a heat source of 100 × 106 m3 at a depth of 5,000 m becomes less than 1 mm/year after 6 years of emplacement. Therefore, the inferred volume of the heat source associated with the 2000 eruption solely corresponds to the magma transported into the shallow depths, while the value given by Miura and Niida (2002) corresponds to the total volume change involved during the 2000 eruption. Terada et al. (2008) also pointed out this inconsistency, and they estimated a heat source volume of about 8 × 106 m3 that is quite similar to our estimate from the observations of heat dis- charge rate near the 2000 eruption vents. Although our estimates of intruded volumes for the 1977 and 1943 eruption sites are comparable to the previous estimates, we cannot rule out that there were no deep intrusions during these two unrests because of the limited geodetic measurements at the time. The discrepancies in parameters of our estimates from other geophysical studies may partly result from the limitations of our thermoelastic model. First, mechanisms other than thermoelastic deflation may also con- tribute to the observed subsidence, although they are not the controlling factor as discussed in section 5. Figure 15. Same as Figure 13 but with the 1943 eruption site. The gray triangles and filled cycles denote the leveling and Interferometric Synthetic Aperture Radar measurements, respectively. 10.1029/2018JB016729 Journal of Geophysical Research: Solid Earth 10.1029/2018JB016729 Journal of Geophysical Research: Solid Earth Fournier & Chardot (JGR, 2012) Wang & Aoki (JGR-Solid Earth, 2019)
34. Monitoring non-magmatic deformation with SAR Ebmeier et al. (J. Appl. Volcanol., 2018) proportion of non-magmatic processes, time as possible (Tilling, 2008). f InSAR and ground based measurement in terms of inferred origins: i) attributed to magmatic processes du d to magmatic processes not associated with eruption, iii) either magmatic or hydrothermal or both in comb tion, iv) attributed to hydrothermal system, v) settling of recent flow deposits and f) displacements associate apse on any scale ü InSAR is more suitable to measure non-magmatic deformation than ground-based measurements. ü Ground-based measurements may not be cost-effective to monitor non-magmatic volcanoes.
35. Hakone Volcano as an analogy to Tatun Volcano? vent area. From the tracks, the average droplet diameter is estimated to have been approximately 2 mm. Since the fall velocity of a droplet of this size is estimated to be approximately 7 m/s (Beard 1977) and the angle between the ash fall tracks and the vertical line was approximately 60°, the wind velocity is estimated to have been 12 m/s (= 7 × tan 60°). We observed no such wind gusts during our stay in the area from 1 h later than this time. Based on this observation, we suspect that the vent proximal area may have been covered by a very dilute base surge, which contained droplets of mud. is is also implied by installed by JMA at Miyagino, in June 30, are the exception and a p 2.8 km above sea level (1.8 km ab nized (Fig. 10c). Mud flow and debris flow from the ve Soon after the recognition of ash Springs Research Institute of (HSRI) formed a joint team with tigate the eruption center and d going down the main stream of Fig. 8 Development of steaming activity at Owakudani before the 2015 Hakone eruption. See Fig. 2 for the shooting lo steaming from the no. 39 SPW (16:48, May 2, 2015; from location c. b Blowout of no. 39 (7:36, May 3, 2015; from near REG around the no. 39 SPW (10:35, May 18, 2015; from near REGMOS). d Ground steaming around the no. 39 SPW (9:04, June camera). Note that no. 39 SPW is not steaming in late June. All photos were taken by HSRI Mannen et al. (Earth Planet. Space, 2018) ü ~70 km from Tokyo, <1 km from the closest resident ü Popular tourist destination ü Well developed hydrothermal systems ü Well monitored ü Long dormance (>600 yrs) before the 2015 eruption Page 3 of 26 nnen et al. Earth, Planets and Space (2018) 70:68 g. 1 Index map of Hakone volcano and Owakudani, center of the 2015 eruption. a Location of Hakone volcano. Yellow triangle indicates active
36. Chronology of the 2001 seismic swarm and 2015 eruption 1 10 100 1000 0 20 40 60 80 0 5 10 15 100 0 100 200 300 a c Baseline Length (cm) Blowout of No.39 SPW in 2015 (May 3) Jan 16, 2015 Mar 4, 2001 Apr 26, 2015 Jul 12, 2001 Nov 12, 2015 Dec 29, 2001 Aug 4, 2015 Sep 20, 2001 Feb 20, 2016 Aug 8, 2002 2015 2001 The 2015 eruption (Jun 29) Blowout of No.52 SPW in 2001(Jul 21) b 2015 2015 2001 2001 d Daily Eq. Number (M ≥ 0) Cumulative DLF Num. ü Similar deformation history ü Similar evolution of earthquake counts ü Different evolution of deep low frequency earthquakes. Mannen et al. (Earth Planet. Space, 2018)
37. Chronology of the 2015 eruption Mannen et al. (Earth Planet. Space, 2018) Page 8 of 26 nen et al. Earth, Planets and Space (2018) 70:68 −6 −5 −4 −3 −2 −1 0 1 2 3 4 5 6 Y, k m 0 1 2 3 4 5 6 7 Depth [km] 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 Depht [km] −2.0 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 2.0 −0.5 0.0 0.5 km km 0 10 20 30 40 50 60 70 80 90100 Time from 22 Apr. [Days] b b c c d d 138.96˚E 139.00˚ 139.04˚ 35.20˚ 35.24˚ 35.28˚N a a KZY KZY KZR KZR C aldera Rim L. Ashinoko L. Ashinoko 15-1 Crater 15-1 Crater Eruption Seismicity Max. g. 6 a Epicenter distribution during the 2015 unrest and eruption of Hakone volcano (from April 22 to July 31). b Depth distribution of a. Seismic- seems to concentrate along planes likely formed by open cracks; an example in the Kojiri area is shown by depth distributions along strike (c) Eruption BO DLF −5 0 5 10 15 Tilt (NS, KZY) −5 0 5 10 15 Tilt (EW, KZY) Seismicity (KZY) 0 20 40 60 80 100 01 11 21 01 11 21 01 11 21 01 11 21 31 0 1000 2000 3000 4000 5000 Cumulative Number April May June July −5 0 5 10 15 Tilt (NS, KZR) −15 −10 −5 0 5 Tilt (EW, KZR) Seismicity (KZR) 0 20 40 60 80 100 01 11 21 01 11 21 01 11 21 01 11 21 31 0 1000 2000 3000 4000 5000 µrad µrad N/hr µrad µrad N/hr Cumulative Number N down E down Fig. 7 Tilt changes (NS and EW components) and seismicity from tiltmeters in Kozukayama (KZY) and Kojiri (KZR). Timings of the initial swarm of the deep low-frequency event (DLF), blowout of the No. 39 SPW (BO) and the eruption are also shown. Seismicity here indicates hourly number of earthquakes within circles shown in Fig. 6. In KZY, seismicity and tilt change are evident in late-April to mid-May and on eruption days. In KZR, both
38. Deformation associated with the 2015 eruption Doke et al. (EPS, 2018) in the Owakudani area were detected by ground-based SAR observation almost at the same time (Doke et al. 2015b). Mannen et al. (2015a) reported that volcanic ashfall started at around 0:30 p.m. (JST) on June 29. e observes surface displacements by multi-mode, right- and left-looking, multi-off nadir angles with a 14-day revisit interval. ALOS-2/PALSAR-2 observed the Hakone Volcano under five different orbits during the period of Fig. 3 Temporal change of local swelling at Owakudani (modified from Doke et al. 2015a). Interferograms are obtained from InSAR pairs in the same orbit with Pair I (Path 18). Color shows wrapped phase difference between 0 and +2π radian. Phase difference of 2π radian corresponds to a slant-range change of 11.9 cm. Contour intervals are 25 m in height Table 3 Open crack and sill (deflation source at depth) parameters estimated from the second inversion a Altitudes are height above sea level of the crack top edges, with ellipsoidal heights in parentheses b East and North coordinates are in WGS84 and refer to the upper left corner of cracks c Parameter was fixed Length (m) Width (m) Altitudea (m) Strike (°) Dip (°) Eastb (°) Northb (°) Opening (cm) 1192.1 298.6 827.6 (868.1) 323.83 88.52 139.030707 35.233585 31.86 261.4 304.1 224.9 (265.4) 0.00 0c 139.027414 35.235381 – 62.10 Fig. 7 Result of the inversion estimating the open crack and deflation sill model. a Modeled surface displacements, and b residuals between model and observation for Pair I. c Modeled surface displacements, and d residuals between model and observation for Pair II. B and C indicate displace- ment areas mentioned in the text. KM and SU show the summits of Kamiyama and Sounzan, respectively. Offset parameters were 1.95 mm for Pair I and –9.20 mm for Pair II. The red line and blue rectangle show the estimated open crack and deflation sill, the parameters of which are listed in Table 3 The observation is modeled by a dike intrusion. Page 12 of 14 ke et al. Earth, Planets and Space (2018) 70:63 ig. 10 Conceptual diagram of the 2015 volcanic activity of Hakone Volcano inferred from the data modeling of InSAR images. a The reservoir of
39. Summary ü Development of volcano geodesy is driven by new observational techniques and new modeling methods. ü Conventional observations and simple and classic modeling are sometimes still useful in understanding the magmatism. ü Sophisticated numerical techniques (e.g., FEM, BEM) are powerful but one needs to understand their limitation as well. ü Deformation signals associated with phreatic eruption are always complicated and resist a simple interpretation. ü Recent unrest of Hakone Volcano could offer insights into the evaluation of future activity of Tatun Volcano.