Articles | Volume 11, issue 1
Research article
22 Mar 2022
Research article |  | 22 Mar 2022

A muographic study of a scoria cone from 11 directions using nuclear emulsion cloud chambers

Seigo Miyamoto, Shogo Nagahara, Kunihiro Morishima, Toshiyuki Nakano, Masato Koyama, and Yusuke Suzuki

One of the key challenges for muographic studies is to reveal the detailed 3D density structure of a volcano by increasing the number of observation directions. 3D density imaging by multi-directional muography requires that the individual differences in the performance of the installed muon detectors are small and that the results from each detector can be derived without any bias in the data analysis. Here we describe a pilot muographic study of the Izu–Omuroyama scoria cone in Shizuoka Prefecture, Japan, from 11 directions, using a new nuclear emulsion detector design optimized for quick installation in the field. We describe the details of the data analysis and present a validation of the results.

The Izu–Omuroyama scoria cone is an ideal target for the first multi-directional muographic study, given its expected internal density structure and the topography around the cone. We optimized the design of the nuclear emulsion detector for rapid installation at multiple observation sites in the field, and installed these at 11 sites around the volcano. The images in the developed emulsion films were digitized into segmented tracks with a high-speed automated readout system. The muon tracks in each emulsion detector were then reconstructed. After the track selection, including straightness filtering, the detection efficiency of the muons was estimated. Finally, the density distributions in 2D angular space were derived for each observation site by using a muon flux and attenuation models.

The observed muon flux was compared with the expected value in the free sky, and is 88 % ± 4 % in the forward direction and 92 % ± 2 % in the backward direction. The density values were validated by comparison with the values obtained from gravity measurements, and are broadly consistent, except for one site. The excess density at this one site may indicate that the density inside the cone is non-axisymmetric, which is consistent with a previous geological study.

1 Introduction

Scoria or cinder cones are a simple volcanic structure, along with stratovolcanoes, shield volcanoes, and lava domes. Understanding the internal structure of scoria cones is important for volcanic hazard assessments. The internal structure of scoria cones has been mainly investigated by geological approaches. Kereszturi and Németh (2012) presented a schematic cross-section of typical scoria cones, and Geshi and Neri (2014) presented detailed photographs of the feeder dike and interior of a scoria cone formed by the 1809 Etna eruption. Yamamoto (2003) investigated outcrops of the interior of scoria cones in the Ojika-jima monogenetic volcano group, Nagasaki Prefecture, Japan. Yamamoto (2003) classified 40 scoria cones according to their degree of interior welding and proposed a link between lava outflow and cone collapse. However, scoria cones with such outcrops are rare, and the internal structure can vary markedly among cones. Therefore, non-destructive methods are required to investigate scoria cones that lack outcrops.

Muography is a non-destructive technique for investigating the internal density structure of large objects, employing the strong penetrating force of muons, which are high-energy elementary particles contained in cosmic rays. Muography has also been used for studying volcanoes, including visualization of a shallow conduit (e.g., Tanaka et al., 2009), detection of temporal changes in water level due to hydrothermal activity (Jourde et al., 2016), and 3D density imaging of a lava dome using a joint inversion of muographic and gravity data (Nishiyama et al., 2017).

In unidirectional muography, the only measurable quantity is the density length, which is the integral of density and length along the muon direction. It has no spatial resolution along the muon path. Therefore, even if an interesting density contrast is found below the crater, this could reflect contributions from other parts of the volcanic body. Similar to X-ray computed tomography, which has been developed as a 3D density imaging technique, muography can obtain 3D spatial resolution by increasing the number of observation directions. In previous studies, muography of volcanoes has been conducted in two or three directions (Tanaka et al., 2010; Rosas-Carbajal et al., 2017). However, the spatial resolution is not sufficient to determine the detailed structure of the volcanic interior. Nagahara and Miyamoto (2018) undertook a 3D density reconstruction based on multi-directional muography and the filtered back-projection technique. Their study showed that it is necessary to increase the number of directions to obtain 3D spatial resolution in volcanological studies.

Nuclear emulsion is a type of muon detector, and has been used for studies of volcanoes (Tanaka et al., 2007; Nishiyama et al., 2014; Tioukov et al., 2019). The trajectories of high-energy charged particles that pass through an emulsion film are recorded as aligned silver grains with micron-scale resolution (Nakamura et al., 2005; Tioukov et al., 2019; Nishio et al., 2020). The positions and slopes of aligned grains in a developed emulsion film are digitized with an automated emulsion readout system (Kreslo et al., 2008; Morishima and Nakano, 2010; Bozza et al., 2012; Yoshimoto et al., 2017). Unlike hodoscopes using scintillator bars (e.g., Saracino et al., 2017) or multi-wire proportional chambers (Oláh et al., 2018), a nuclear emulsion film does not have temporal resolution. In contrast, an emulsion detector does not require electricity, which facilitates the installation of such detectors around volcanoes where infrastructure is not well developed.

In muographic studies of a volcano, contamination by low-momentum particles must be removed to derive the correct density (Nishiyama et al., 2014, 2016). Thus, nuclear emulsion detectors have often been used as an emulsion cloud chamber (ECC), which comprises alternating layers of films and lead or iron plates (e.g., Kodama et al., 2003). An ECC detector can measure the momentum of the charged particles, one by one, by detecting deflection angles caused by multiple Coulomb scattering (Agafonova et al., 2012). For multiple Coulomb scattering, there is a relationship between the maximum detectable momentum pmax and position resolution yreso as follows using the first term of Eqs. (33.16) and (33.20) in Tanabashi et al. (2018):

(1) p max α X 0 - 0.5 x 1.5 y reso - 1 ,

where α is a constant, X0 is the radiation length of a material, and x is the thickness of the material. The position resolution of the newest scintillator hodoscope or MWPC is on the order of 1 mm (Saracino et al., 2017; Oláh et al., 2018). In the case of nuclear emulsion, the resolution is about 1 µm. When using ECC, the thickness of the material can be reduced to 1/100 while maintaining the same pmax, which is advantageous in terms of transportation in the field.

A new design of the ECC detector was also required for its rapid installation at multiple observation sites in the field. In a previous study of volcano observations using the ECC detector (Nishiyama et al., 2014), rapid installation of the detector was not required because the number of observation sites was just one. It is also important to establish a data analysis procedure for the muon tracks recorded by the ECC detectors. To derive an accurate density value for the volcanic body, it is necessary to remove low-momentum contamination, estimate the detection efficiency, and validate the results. In addition, for bias-free 3D imaging by multi-directional muography, the installed muon detectors must show similar performance.

2 Izu–Omuroyama scoria cone

The Izu–Omuroyama scoria cone (345411′′ N, 1390540′′ E; 580 m a.s.l.) is one of the largest scoria cones in the world, and is part of the Higashi Izu monogenetic volcano group (Aramaki and Hamuro, 1977), which is located in the northeastern Izu Peninsula, Ito City, Shizuoka Prefecture, Japan. It is considered to have formed about 4000 years ago, based on 14C dating (Saito et al., 2003). The basal diameter is 1000 m, the height is 280 m from the base, and the typical slope of its flanks are 29–32. The center of the cone contains a crater that is 250 m wide and 40 m deep. The volume of the cone is 71 × 106 m3, and lava with a volume of ∼108 m3 has flowed out from the base of the cone (Koyano et al., 1996). The lava is a basaltic andesite of 54 wt %–56 wt % SiO2 (Hamuro, 1985).

Although the shape of the Izu–Omuroyama scoria cone appears to be axisymmetric (Fig. 1), a geological study suggested it has an anisotropic structure due to the following reasons. (i) During/after the growth of the cone, some interior parts became welded due to loading, residual heat, and a low cooling rate. As a result, some denser material formed. (ii) At the end of the eruption, a lava lake was formed in the crater, and the lava flowed out to the western foot of the cone. (iii) There is a small crater on the south side of the cone, which is thought to have formed when the main crater was blocked at the end of the eruption (Koyano et al., 1996).

Figure 1Photograph of the Izu–Omuroyama scoria cone from the northwest, taken by an unmanned aerial vehicle (Koyama, 2017).


The bulk density of typical continental crust is about 2.6–2.7×103 kg m−3. The bulk densities reported for scoria deposits are 0.84–1.01×103 kg m−3 (Taha and Mohamed, 2013) and 0.56–1.20×103 kg m−3 (Bush, 2001). Therefore, the maximum expected density contrast is about 1.4–2.0×103 kg m−3, due to the difference in porosity between welded rocks and scoria deposits. In addition, the Izu–Omuroyama scoria cone is an ideal target for multi-directional muography due to the accessibility of detector sites and absence of muographic shadows from any direction caused by other topographic features.

3 Multi-directional muography observations using emulsion cloud chambers

3.1 Detector design

Emulsion films were manufactured by pouring 70 µm of nuclear emulsion on both sides of a 180 µm thick plastic base. The size of a film is 125 × 100 mm. The films were vacuum-packed in a light-blocking envelope to maintain their planar form, which prevented air bubbles from forming between the envelope and film, and made it easy to handle the films in the field.

The detector used for the 2018 observations is basically the same as that of Nishiyama et al. (2014), and only the number of lead plates was different. The former consists of 20 films and 9 plates of 1 mm thick lead, the latter consists of 20 films and 19 lead plates. At the time of installation in 2018, the films, lead plates, and supports were all in pieces and, therefore, much time and effort was required for assembly in the field. A more efficient detector design was required for rapid and error-free installation.

The detector used in the 2019 observations was improved. It consists of an ECC and an outer box. The ECC consists of 20 emulsion films and 19 lead plates, each 1 mm thick (Fig. 2a). An aluminum frame was fixed to a lead plate with a thin sheet of glue, and then an emulsion film with the light-blocking envelope was attached with scotch tape. In this paper, we term this unit the emulsion–lead plate (EL plate; Fig. 2a). The EL plate was designed for quick assembly in the field.

Figure 2Design of the ECC and outer box. (a) Schematic cross-section of the EL plates and an ECC. The EL plates consist of a 1 mm thick aluminum frame, 1 mm thick lead plate, 100 µm thick glue sheet that fixes a lead plate to an aluminum frame, and an emulsion film with a light-blocking envelope. An ECC consists of 19 EL plates and an emulsion film with a plastic plate. (b) Schematic of the aluminum outer box. The thickness of the aluminum plate is 10 mm. The ECC shown in (a) was set inside this box. There are four holes for feed screws in the front plate.


The outer box consists of 10 mm thick aluminum plates (Fig. 2b). The outer size of this box is 190 mm in width, 155 mm in height, and 90 mm in depth. An ECC and strong springs were placed in the box. There are four screw holes on one side of the box, and by turning the bolts and pushing the spring plate, a uniform pressure ( 105 Pa) was applied to the ECC. This pressure prevents the film from stretching and shrinking due to temperature changes.

Given that there is no temporal resolution in emulsion films, ordinary ECC detectors cannot distinguish whether cosmic ray tracks pass the ECC during muographic observation or during transportation and standby. Thus, we also added a similar feature to previous muographic studies using emulsion films. Researchers have used emulsion films with a different alignment during muon observation and during standby (e.g., Tanaka et al., 2007). In the present study, the corners of the EL plates were aligned during the muon observations, while the corners were intentionally shifted a few millimeters horizontally and fixed with clamps during standby (Fig. 3). This alignment difference distinguishes passing charged particles during non-observation and observation periods by pattern matching each emulsion film. By using this procedure, the time to set the alignment between each EL plate in the field is <30 s. Although the muon tracks that pass through an ECC during the alignment set-up may become noise, our procedure reduced such tracks.

Figure 3(a) View of the EL plates from above during standby. The EL plates were intentionally shifted a few millimeters horizontally and fixed with a pair of steel plates and a clamp. The red lines represent the muon tracks in this alignment. (b) View from above during the observations. The EL plates were aligned to the side of the outer box, and fixed by the springs and feed screws. The blue lines represent muon tracks during observations. Note that the red tracks cannot be reconstructed in this alignment.


3.2 Installation

The detectors were installed at three sites in 2018 and eight sites in 2019 around the Izu–Omuroyama scoria cone (Fig. 4; Table 1). Each detector was buried in a hole that was about 40 cm deep to avoid high temperatures due to direct sunlight. This is done because the number of latent image specks decreases and the number of randomly generated specks increases under high-temperature conditions (Nishio et al., 2020).

Figure 4Topography of the Izu–Omuroyama scoria cone. White dots represent observation sites.


Table 1Effective area and muon exposure time for each detector.

Download Print Version | Download XLSX

The installation procedure at each observation site in 2019 was as follows (Fig. 5).

  1. Carry the outer box and EL plates to the observation site.

  2. Measure the coordinates of the site with a hand-held GPS (Garmin; model GPS eTrex 30J). The typical uncertainty of the latitudinal and longitudinal coordinates is 3 m.

  3. Dig a hole in the ground with horizontal dimensions of 60 × 40 cm and a depth of 40 cm.

  4. Flatten the base of the hole, place a plastic bag inside the hole, and lay down a piece of plywood.

  5. Put double-sided tape on the bottom of the outer box and place it on the plywood.

  6. Put the stack of EL plates into the box and quickly align them (<30 s).

  7. Close the cap of the outer box.

  8. Turn the feed screws to increase the pressure.

  9. Measure the attitude of the outer box (i.e., the yaw, that is, absolute azimuth angle; roll; and pitch). The yaw is measured with a fiber optic gyro (Japan Aviation Electronics Industry Ltd.; model FOG JM7711; Watanabe et al., 2000), and the roll and pitch are measured by the digital leveler. The typical errors in the yaw, roll, and pitch are 8.7×10-3, 1.0×10-3, and 1.0×10-3 radians, respectively.

  10. Cover with styrofoam to avoid heating from the ground surface.

  11. Close the plastic bag to keep water out.

  12. Backfill the hole.

The time taken for this installation was  2 h for each site, and we installed detectors as three sites in a day in 2019. The detector retrieval procedure was the opposite of the installation procedure. The 380 films were developed in a darkroom. The deposited silver particles on the surface of the films were removed with anhydrous ethanol. The gelatin of the sensitive layer was swollen with a glycerin solution to obtain the optimum thickness for an automated track readout system, which is described in the next section.

Figure 5Photographs showing the installation procedure. (a) Dig a hole and place a plywood sheet in the bottom. (b) Place the outer box in the hole and put a stack of EL plates into the box. The plates were aligned over a period of <30 s. After closing the top plate of the box, the feed screws were tightened to increase the pressure. (c) The yaw, roll, and pitch were measured with a fiber optic gyro (FOG) and digital leveler.


4 Track reconstruction, selection, and detection efficiency estimation

4.1 Track reconstruction

A track of a high-energy charged particle is recorded as an aligned line of silver grains in the emulsion film (e.g., Nakamura et al., 2005). The images in the 380 nuclear emulsion films were scanned and the positions and slopes of the tracks were digitized by the “HTS”, a high-speed automated track readout system at Nagoya University (Yoshimoto et al., 2017). For each ECC, the tracks of the charged particles were digitally reconstructed from the segmented tracks in 20 films. NETSCAN 2.0 software was used for track reconstruction (Hamada et al., 2012). NETSCAN 2.0 rapidly corrects for film distortions and local misalignments between films by using many tracks recorded over a large area. It then outputs all possible connections as the final result. NETSCAN 2.0 has been used in various fields, such as neutrino physics (Hiramoto et al., 2020), cosmic ray astronomy (Takahashi et al., 2015), and muographic studies of Egyptian pyramids (Morishima et al., 2017). The typical procedure for the track reconstruction is as follows.

  1. Reconstruct the “base track”, which is connected between the emulsion layers across the plastic base of 170 µm in a film.

  2. Reconstruction of the “linklet”, which is the base track pair between adjacent films across lead plates.

  3. Reconstruction of the tracks that connect across the whole ECC. If no base track was found in two consecutive films on the extension of a track, then the track was considered to have stopped.

For example, in ECC_ID = 02, 8.9 × 106 base tracks, 3.2 × 106 linklets in a pair of adjacent films, and 1.7 × 107 tracks in an entire ECC were reconstructed.

4.2 Track selection

NETSCAN 2.0 outputs all possible track connections. Therefore, it is necessary to carefully select the tracks for the muographic analysis. A schematic example of the output tracks is shown in Fig. 6. Most of the branches can be considered to represent contamination by fake base tracks caused by random noise, or the coincidental occurrence of low-energy positrons/electrons on parallel slopes in the vicinity of the real tracks (e.g., Fig. 6; cases 2 and 3). For reference, the position dependence of noise density is described in Appendix A. Some branches consist of a pair of straight tracks with small closest distances and similar angles (Fig. 6; case 4). In this case, the two tracks should be separated.

Figure 6Schematic examples of typical reconstructed tracks in an ECC obtained by NETSCAN 2.0. Upstream means towards the volcanic cone side and downstream means the backward free sky direction. Case 1: a straight track without any branches. Case 2: a straight track with a branch in the middle and downstream films. The track branch in the middle was rejected by selection step 1. The branch in the most downstream film was merged into the straight track by selection step 3a. Case 3: branches in the upstream and middle films. Both branches were merged into a straight track by selection step 3a. Case 4: a pair of straight tracks with small closest distances and similar angles. If the shared proportion of the track length was <20 %, the tracks were divided into two different tracks by selection step 3b.


The following χ2/ ndf value was calculated for all tracks for the low-momentum cut-off:

(2) χ 2 / ndf = m Δ θ R m σ R m 2 + Δ θ L m σ L m 2 / ndf ,

where ndf is the number of degrees of freedom and m is the index of adjacent film pairs (i.e., [1,2], [2,3], [3,4], …, and [18,19], [19,20] in Fig. 6) or with one skip if there was a base track inefficiency (i.e., [1,3], [2,4], [3,5], …, [17,19], [18,20]). ΔθRm=(Δθxm×tanθx+Δθym×tanθy)/tan2θx+tan2θy and ΔθLm=(Δθym×tanθx-Δθxm×tanθy)/tan2θx+tan2θy, and Δθxm and Δθym are angular differences along the x, y coordinates of the ECC. σRm and σLm are the root mean square of ΔθRm and ΔθLm, calculated for every adjacent film pair in every ECC (Fig. 7). Figure 8 shows the distribution of χ2/ ndf for all tracks in an ECC.

Figure 7Examples of σR and σL as a function of tan2θx+tan2θy. The values were determined by the ECC and used to calculate the value of Eq. (2).


Figure 8Example of the χ2/ ndf distribution for selected tracks as a function of θ=tan2θx+tan2θy in an ECC. (a) 0θ<0.4, (b) 0.4θ<0.8, (c) 0.8θ<1.2, and (d) 1.2θ<1.6.


The procedure for track selection is as follows.

  • 1.

    Select tracks that start from one of the two most upstream (i.e., summit cone side) films and stop at one of the two most downstream films.

  • 2.

    Select tracks with χ2/ndf<5.0.

  • 3.

    If a track has any branches, then do the following:

    • a.

      If the shared proportion of track length is  20 %, choose the longest branch. If the track lengths are the same, then choose the branch with the smallest χ2/ ndf.

    • b.

      If the shared proportion of track length is <20 % (Fig. 6; case 4), then divide the branches into two tracks.

We estimated the effect of the straightness filtering using χ2/ndf<5.0. Figure 9 shows the momentum filtering efficiency. This figure was derived from a simple simulation in which the interaction of charged particles inside the ECC was assumed to be multiple Coulomb scattering only, and the scattering angle was approximated by a Gaussian distribution. The path length in the lead plates (i.e., x in Eq. 1) becomes longer when the track has a larger slope, and thus the maximum detectable momentum (i.e., pmax in Eq. 1) also becomes higher. Based on the background noise study by Nishiyama et al. (2016), the size of the mountain body used in the simulation and the Izu–Omuroyama scoria cone is broadly the same, and thus the rejection efficiency should be sufficient. For example, after the track selection, 1.7 × 106 tracks were selected at the site N.

Figure 9Survival ratio of muons after the straightness cut-off as a function of momentum. (a) Track angles with tan (relative azimuth)=0.0 and tan (elevation angle)=0.0. (b) Track angles with tan (relative azimuth)=1.0 and tan (elevation angle)=1.0. The path length in the lead plates becomes longer when the track has a larger slope, and thus the remaining momentum also becomes higher for the latter case. The momentum values at a survival ratio of 0.5 are 0.6 and 1.1 GeV / c, respectively.


4.3 Detection efficiency estimation

The muon detection efficiency can be estimated by investigating the percentage of tracks that have a base track in a film. In this paper, we term this percentage the “fill factor”. The fill factor ε can be defined as follows:

(3) ε j θ x , θ y = N j θ x , θ y N j - 1 , j + 1 θ x , θ y ,

where j is a film ID, Nj-1,j+1θx,θy is the number of tracks in which base tracks were found in films j−1 and j+1, and Nj(θx, θy) is the number of tracks in which base tracks were found in films j−1, j, and j+1. The fill factor depends on the films and track slopes θx and θy. The position dependence of the fill factor is described in Appendix A.

Using the fill factor εj(θx,θy) and εjθx,θy=1-εjθx,θy, the muon detection efficiency ϵ in an ECC can be calculated as follows:

(4) ϵ θ x , θ y = hit pattern ε 1 × ε 2 × ε 3 × × ε 18 × ε 19 × ε 20 ,

where “hit pattern” is the summation for all possible hit patterns (e.g., ε1×ε2×ε3××ε18×ε19×ε20 or ε1×ε2×ε3××ε18×ε19×ε20) from the track selection conditions described in section 4.2 (Fig. 10). An example of the angular distribution of the fill factor εj(θx, θy) and muon detection efficiency ϵ(θx, θy) in an ECC is shown in Fig. 11. The statistics of observed muons were limited in some angular bins by the thick volcanic cone. However, the statistics were sufficient in the backward region (i.e., elevation angle <0.0). We used the distribution of the negative elevation angular region instead of the positive region, because it has enough statistics and the optics of the HTS has an approximately two-fold rotational symmetry.

Figure 10Example of all hit patterns and the products of fill factors in Eq. (4) when base tracks are found in film ID numbers 4 and 8. The red lines indicate the reconstructed tracks and the short black lines represent the base tracks found in the films.


Figure 11Examples of the angular distribution of the fill factor in some films and the efficiency of an ECC. (a) Fill factor for film ID = PL01 (most upstream film) and ECC ID = 02 at site N. (b) Fill factor for film ID = PL02 and ECC ID = 02. (c) Muon detection efficiency for ECC ID = 02 as evaluated by Eq. (4). The horizontal axis is the tangent of the azimuth angle; the vertical axis is the tangent of the elevation angle; the colors represent the fill factor/efficiency values. A positive elevation angle means the muon path is from the cone; a negative elevation angle means the muon path is from the backward free sky. The gray color means there were no observed muons in the angular bin due to the thick volcanic cone.


5 Results

The average density along muon path was determined for each observation site. We used the muon flux model of Honda et al. (2004), energy loss model of Groom et al. (2001), and topography around the Izu–Omuroyama scoria cone from the Geospatial Information Authority of Japan (, last access: 11 March 2022). The coordinates of the observation site, direction, sensitive area, thickness of the ECC detectors, and observation time were used to calculate the expected number of muons at each observation site. The expected number of muons can be calculated as a function of the average density ρk along the path:

(5) N k simu ( ρ k ) = f k ρ k L k × S k × Ω k × T × ϵ k ,

where k is the index of an angular bin, fk(ρkLk) is the penetrating muon flux (calculated from the muon flux model, energy loss model, and path length Lk), Sk is the sensitive area of the ECC, Ωk is the solid angle, ϵk is the muon detection efficiency, and T is the observation time.

The angular bin size used for calculating the expected value was (0.01)2 in terms of the tangent. The angular bins were then merged to improve the statistical accuracy of the observed values. This merging procedure is useful in topology where a small change in elevation angle can dramatically change the path length in the volcano. If k is the index of the angular bins of (0.01)2 and the bins belong to a larger angular bin i, then the following equation holds:

(6) N i merged ρ i = k N k simu ρ i ,

where ρi is the density of the merged angular bin i. If Niobs is the number of the detected muons in the angular bin i, then we can uniquely determine the density value ρi, such that Nimergedρi=Niobs. The lower limit ρilow and upper limit ρiup caused by the statistical error in Niobs can also be estimated as follows:


An example of the derived density map is shown in Fig. 12. All results are shown in Figs. B1–B10 (Appendix B).

Figure 12Data for observation site N. (a) Map, topography, and viewing range; (b) path length of the volcanic cone; (c) muon counts Niobs; (d) density ρi. The maximum value of the color bar indicates a density of >2.5×103 kg m−3 and the minimum value is <0.5×103 kg m−3(e) Density error Δρ=ρiup-ρilow/2. The maximum value of the color bar indicates a density error of >0.25×103 kg m−3.


The definition of the angular bin areas was based on the following. The size of the angular bins was (0.2)2 when the elevation angle was 0.1 to 0.5 in tangent terms. When the elevation angle was >0.5, the angular bin size was (0.1)2. If the observed muon count in the bin was <25, the angular bin was manually merged with adjacent bins to improve the statistical error. The angular bin with a near-surface path (path length < 30 m) was excluded to avoid ambiguity between the actual topography and digital elevation map. The attitude errors of each muon detector also contribute to the path length ambiguity, especially near the surface of the cone.

6 Validation

Firstly, we validated the observed muon flux by comparing it with the muon flux model in the free sky region. The average and standard deviation of the ratio between the sites were 88 % ± 4 % in the forward direction and 92 % ± 2 % in the backward direction, except for the NNE site (Fig. 13). An example of observed/expected muon flux ratio angular distribution of the site N is shown in Fig. 14. As can be seen in this figure, in each detector site, a non-uniform distribution of the observed/expected muon flux ratio exists. Such deviations were 4 %–7 % except for the forward directions at the sites SE and NNE. For reference, a 10 % error in the flux corresponds to a 4 % error in the density length at a tan(elevation angle) = 0.2 and density length = 1000 m (water equivalent). These deviations were less than the errors caused by the muon statistics. The discrepancy for the NNE site is discusses in the next section.

Figure 13The observed / expected muon flux ratio for each observation site in the free sky region. The plot represents the average value of the ratio in tangential angular space, and the error bars are the standard deviations at each site.


Figure 14Example of the observed / expected muon flux ratio in the free sky region at site N. The horizontal axis represents the azimuth angle, and the vertical axis represents the elevation angle. Positive elevation angle means the muons come from forward directions (the volcanic cone). Negative elevation angle means the muons come from backward free sky directions. The typical deviation of the ratio is 4 %–7 % in each site.


Secondly, we compared the density of the entire volcanic cone determined by gravity data with that obtained by muography. Table 2 shows the density determined from each observation site when the cone is assumed to be uniform. The calculation of the overall density ρ is as follows:

(9) ρ = i ρ i V i i V i ,

where i is the index of the angular bins and Vi is the volume of the volcanic cone cut off by the angle bin i. Based on the gravity study of Nishiyama et al. (2020), the density of the Izu–Omuroyama scoria cone is 1.39±0.07×103 kg m−3. The overall density derived by muography at each observation site is 1.42–1.53 × 103 kg m−3, except for one site. These values are broadly consistent with the density determined from gravity data, except for the observation site W (1.72 × 103 kg m−3).

Table 2Overall bulk density obtained by muography, assuming that the density is uniform in the volcanic cone.

Download Print Version | Download XLSX

7 Discussion

For the observed / expected muon flux ratio in the free sky region, the values in the forward direction are less than in the backward direction at many observation sites. This could be because the detectors were buried in holes on steep slopes ( 30), and our analysis might not account for that effect. Due to the steep slope, muons arriving from the forward direction need to penetrate some amount of soil, whereas muons from the backward direction can enter the detector without being affected by the soil cover. In addition, the resolution of the detector coordinates is  3 m, which might also contribute to the discrepancy.

Some density results from near the ground surface are complex. Some regions near the path length of 30 m appear to have relatively higher or lower density than the other data (e.g., Figs. B6, B9). One possible reason for this is the error in the detector attitude. Near the surface of the volcanic cone, the difference between the calculated and actual path lengths may become larger due to the error in the detector attitude.

The anomalous data for the NNE site also warrants further consideration. The reason for this might be a difference between the digital elevation map and actual topography. There is a stone wall in front of the buried detector at this site, which is about 1 m high and located on the volcanic cone side. The grid size of the digital elevation map used in this study is 1 m, and thus the map might not record this steep gradient.

In summary, errors in the position and attitude of the detectors, and the accuracy of the DEM, might cause a misfit between the DEM and actual topography. These are the main reasons for the discrepancy between the observed and expected muon flux.

The discrepancy between the observed and expected muon flux was ±4 % in the forward direction and ±2 % in the backward direction between the detectors. In addition, the typical deviation inside each site was 4 %–7 %. These values are smaller than the statistical error of the observed muons used to determine the density of the volcanic cone, and thus they were not significant for our observations. It is interesting to consider if an improvement in the accuracy of the detector position and attitude, and the DEM, would decrease this systematic error. For example, the ±4 % deviation in the forward direction would be expected to decrease to ±2 %, because the misfit effect is less in the backward direction. Further improvements will require simulation of the expected muon flux that take into account more processes and verification of the systematic errors associated with the ECC detectors.

The obtained density values (1.42–1.53×103 kg m−3; this study) and 1.39±0.07×103 kg m−3 (Nishiyama et al., 2020) for the Izu–Omuroyama scoria cone are broadly consistent (Table 2). In a study by Rosas-Carbajal et al. (2017), the density value obtained from gravity data was 0.5 × 103 kg m−3 larger the value obtained by muography. In our validation, this discrepancy does not exist. As Rosas-Carbajal et al. (2017) suggested, the discrepancy might be due to differences in the filtering performance for low-momentum particles, as shown in Fig. 9.

The higher density obtained at site W cannot be explained by the systematic errors described above. One possible reason for this is an actual high-density structure in front of the site. This hypothesis is consistent with the fact that lava flowed out from the crater lake to the west (Koyano et al., 1996).

8 Conclusions

A muographic study of the Izu–Omuroyama scoria cone was undertaken in 11 directions. The ECC detector design was optimized for quick installation in the field. We mounted the 11 detectors beneath the ground, surrounding the volcanic cone. The tracks of charged particles that passed through the ECCs were reconstructed using the automated emulsion track readout system HTS and NETSCAN 2.0 software. After track selection, including momentum filtering and efficiency estimation, the density profiles in 2D angular space were derived for each observation site. The methods described in this paper can be applied to the observation of other volcanoes and target objects.

We compared the observed muon flux to the expected value from a muon flux model in the free sky region. The muon flux difference between each detector was 4 % in the forward directions and 2 % in the backward directions, and the typical deviations in each site were about 4 %–7 %. The errors in the detector coordinates and attitude, and DEM, are the main cause of the discrepancy between the observed and expected muon flux.

In addition, we also compared our results with the overall volcanic cone density estimated from gravity data, which are broadly consistent, apart from the W site. This discrepancy for the W site cannot be explained by the systematic errors discussed in the previous section and the statistical error of the observed muons. It might reflect a high-density structure located in the western flank of the volcano. Further 3D density reconstructions of the Izu–Omuroyama scoria cone are ongoing using the data set described in this paper.

Appendix A

We here consider how the position dependency of the detected tracks affects the density results.

The fill factor of the base tracks (i.e., track detection efficiency in a film) also depends on the position of the scanned film. The typical causes of the decrease of the fill factor are heterogeneous thickness of the emulsion layers, some dusts or scratches on the emulsion surface, and the poorly tuned parameters for the scanning.

Figure A1 shows the position distribution of the fill factor of all films of an ECC. For example, on the upper right, the films tend to have a low fill factor (e.g., a–f, h, k, l, q). This part has a larger thickness of the emulsion layer because the drops of the glycerin solution were left in the upper right corner when drying after soaking due to the structure of the developing racks. Figure A1s and t show a larger area of the low fill factor on the right and left. The reason might be the poorly tuned parameters used for the scanning.

Compared to the size of the scoria cone, the ECC is a very small “element”; thus the local position dependence of the fill factor can be approximately treated as an average fill factor εj(θx,θy). The decrease of the fill factor is reflected in εj(θx,θy) in Eq. (4). Finally, εj(θx,θy), which encompasses the effects of the local decrease of the fill factor, is effectively used to derive the angle-dependent muon detection efficiency.

We should also consider the position dependency of noise. Local high density of random silver grains caused by any chemical conditions, or fake images produced by scratches on the films tend to create a group of fake tracks concentrated in one place. In addition, such fake tracks tend to have small slopes by scanning with automated emulsion readout system. If such noise is continuous at the same location on the film, they might make many parallel tracks at a certain slope and give a systematic error in the result. However, such possibility has been eliminated by the track selection algorithm described in the Sect. 4.2, because such concentrated tracks in position and angular space make numerous branches in track connection (see Fig. 6) and tracks with such branches were removed in the selection procedure. Figure A2 shows the number of selected tracks with small slope per mm2 in each observation site. The histogram roughly fits the Poisson distribution and there is no remarkable excess. The difference in the peak in the histograms depends on the difference in the exposure time (SE, W, NNE), existence of topography in the backward direction (NE), pitch angle of the detector attitude (e.g., SW has larger pitch angle, thus less tracks of the small slopes from the backward direction), and the difference in muon detection efficiency.

Figure A1The position distribution of the fill factor in each film of ECC02. (a)(t) represent PL01–PL20, respectively.


Figure A2(a) The position distribution of the number of selected tracks per mm2 in the ECC02. (b–l) The number of selected tracks per mm2 (the black line) with the fitting result of a Poisson distribution (=αμxe-μ/x, the red line) in each observation site, respectively. The tracks selected for this figure came from in the backward direction and have small slopes (|tanθx|<0.5 and |tanθy|<0.5).


Appendix B

The density results for each observation site are shown in Figs. B1–B10.

Figure B1Observation site NW. (a) Map, topography, and viewing range; (b) path length of the volcanic cone; (c) muon counts Niobs; (d) density ρi. The maximum value of the color bar indicates a density of >2.5×103 kg m−3 and the minimum value is <0.5×103 kg m−3. (e) Density error Δρ=ρiup-ρilow/2. The maximum value of the color bar indicates a density error of >0.25×103 kg m−3.


Figure B2Observation site W. (a) Map, topography, and viewing range; (b) path length of the volcanic cone; (c) muon counts Niobs; (d) density ρi. The maximum value of the color bar indicates a density of >2.5×103 kg m−3 and the minimum value is <0.5×103 kg m−3. (e) Density error Δρ=ρiup-ρilow/2. The maximum value of the color bar indicates a density error of >0.25×103 kg m−3.


Figure B3Observation site WSW. (a) Map, topography, and viewing range; (b) path length of the volcanic cone; (c) muon counts Niobs; (d) density ρi. The maximum value of the color bar indicates a density of >2.5×103 kg m−3 and the minimum value is <0.5×103 kg m−3. (e) Density error Δρ=ρiup-ρilow/2. The maximum value of the color bar indicates a density error of >0.25× 103 kg m−3.


Figure B4Observation site SW. (a) Map, topography, and viewing range; (b) path length of the volcanic cone; (c) muon counts Niobs; (d) density ρi. The maximum value of the color bar indicates a density of >2.5×103 kg m−3 and the minimum value is <0.5×103 kg m−3. (e) Density error Δρ=ρiup-ρilow/2. The maximum value of the color bar indicates a density error of >0.25×103 kg m−3.


Figure B5Observation site S. (a) Map, topography, and viewing range; (b) path length of the volcanic cone; (c) muon counts Niobs; (d) density ρi. The maximum value of the color bar indicates a density of >2.5×103 kg m−3 and the minimum value is <0.5×103 kg m−3. (e) Density error Δρ=ρiup-ρilow/2. The maximum value of the color bar indicates a density error of >0.25×103 kg m−3.


Figure B6Observation site SE. (a) Map, topography, and viewing range; (b) path length of the volcanic cone; (c) muon counts Niobs; (d) density ρi. The maximum value of the color bar indicates a density of >2.5×103 kg m−3 and the minimum value is <0.5×103 kg m−3. (e) Density error Δρ=ρiup-ρilow/2. The maximum value of the color bar indicates a density error of >0.25×103 kg m−3.


Figure B7Observation site ESE. (a) Map, topography, and viewing range; (b) path length of the volcanic cone; (c) muon counts Niobs; (d) density ρi. The maximum value of the color bar indicates a density of >2.5×103 kg m−3 and the minimum value is <0.5×103 kg m−3. (e) Density error Δρ=ρiup-ρilow/2. The maximum value of the color bar indicates a density error of >0.25×103 kg m−3.


Figure B8Observation site E. (a) Map, topography, and viewing range; (b) path length of the volcanic cone; (c) muon counts Niobs; (d) density ρi. The maximum value of the color bar indicates a density of >2.5×103 kg m−3 and the minimum value is <0.5×103 kg −3. (e) Density error Δρ=ρiup-ρilow/2. The maximum value of the color bar indicates a density error of >0.25×103 kg m−3.


Figure B9Observation site NE. (a) Map, topography, and viewing range; (b) path length of the volcanic cone; (c) muon counts Niobs; (d) density ρi. The maximum value of the color bar indicates a density of >2.5×103 kg m−3 and the minimum value is <0.5×103 kg m−3. (e) Density error Δρ=ρiup-ρilow/2. The maximum value of the color bar indicates a density error of >0.25×103 kg m−3.


Figure B10Observation site NNE. (a) Map, topography, and viewing range; (b) path length of the volcanic cone; (c) muon counts Niobs; (d) density ρi. The maximum value of the color bar indicates a density of >2.5×103 kg m−3 and the minimum value is <0.5×103 kg m−3. (e) Density error Δρ=ρiup-ρilow/2. The maximum value of the color bar indicates a density error of >0.25×103 kg m−3.


Code availability

The participants of this study did not agree for their code to be shared publicly, so the supporting code is not available.

Data availability

The data of this study are available from the corresponding author upon reasonable request.

Author contributions

SM, SN, MK, and YS discussed the significance of the study and were involved in the study design. KM checked the quality of nuclear emulsion, produced the films, and prepared the developing solution. SM and SN designed the new emulsion film detector, installed the detectors, and developed the films after muon exposure. TN operated the automated track readout system for the scanning of films. SM and SN analyzed the track data and determine the density values. All of authors considered the results of the study.

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The authors thank Hideaki Aoki and his colleagues of Ike-kankou for collaborating on our study. We also thank Masakazu Ichikawa of the Earthquake Research Institute, University of Tokyo, for support during the observation campaign. We are also grateful for the technical support of the staff and students in F-lab, Nagoya University, especially with the nuclear emulsion films. This research was supported by JSPS KAKENHI Grant 19H01988, an Izu Peninsula Geopark Academic Research Grant (2018), the Joint Research Program of the Institute of Materials and Systems for Sustainability at Nagoya University (2017–2021), and a JSPS Fellowship (Nagahara; Grant DC2, 19J13805).

Review statement

This paper was edited by Ralf Srama and reviewed by two anonymous referees.


Agafonova, N., Aleksandrov, A., Altinok, O., et al.: Momentum measurement by the multiple Coulomb scattering method in the OPERA lead-emulsion target, New J. Phys., 14, 013026,, 2012. 

Aramaki, S. and Hamuro, K.: Geology of the Higasi–Izu monogenetic volcano group, Bull. Earthq. Res. Inst. Univ. Tokyo, 52, 235–278, (last access: 11 March 2022), 1977 (in Japanese). 

Bozza, C., D'Ambrosio, N., De Lellis, G., et al.: An integrated system for large scale scanning of nuclear emulsions, Nucl. Instrum. Methods Phys. Res. A, 703, 204–212,, 2012. 

Bush, A. L.: Construction Materials: Lightweight Aggregates, in: Encyclopedia of Materials: Science and Technology, Elsevier, 1550–1558,, 2001. 

Geshi, N. and Neri, M.: Dynamic feeder dyke systems in basaltic volcanoes: the exceptional example of the 1809 Etna eruption (Italy), Front. Earth Sci., 2, 13,, 2014. 

Groom, D. E., Mokhov, N. V., and Striganov, S. I.: Muon stopping power and range tables 10 MeV–100 TeV, At. Data Nucl. Data Tables, 78, 183–356,, 2001. 

Hamada, K., Fukuda, T., Ishiguro, K., Kitagawa, N., Kodama, K., Komatsu, M., Morishima, K., Nakano, T., Nakatsuka, Y., Nonoyama, Y., Sato, O., and Yoshida, J.: Comprehensive track reconstruction tool “NETSCAN 2.0” for the analysis of the OPERA Emulsion Cloud Chamber, J. Instrum., 7, P07001,, 2012. 

Hamuro, K.: Petrology of the Higashi–Izu monogenetic volcano group, Bull. Earthq. Res. Inst. Univ. Tokyo, 60, 335–400, (last access: 11 March 2022), 1985. 

Hiramoto, A., Suzuki, Y., Ali, A., et al.: First measurement of νμ and νμ charged-current inclusive interactions on water using a nuclear emulsion detector, Phys. Rev. D, 102, 072006,, 2020. 

Honda, M., Kajita, T., Kasahara, K., and Midorikawa, S.: New calculation of the atmospheric neutrino flux in a three-dimensional scheme, Phys. Rev. D, 70, 043008,, 2004. 

Jourde, K., Gibert, D., Marteau, J., d'Ars, J. B., and Komorowski, J. C.: Muon dynamic radiography of density changes induced by hydrothermal activity at the La Soufrière of Guadeloupe volcano, Sci. Rep., 6, 33406,, 2016. 

Kereszturi, G. and Németh, K.: Monogenetic Basaltic Volcanoes: Genetic Classification, Growth, Geomorphology and Degradation, Updates in Volcanology – New Advances in Understanding Volcanic Systems, Karoly Nemeth, IntechOpen,, 2012. 

Kodama, K., Hoshino, K., Komatsu, M., Miyanishi, M., Nakamura, M., Nakamura, T., Nakano, T., Narita, K., Niwa, K., Nonaka, N., Sato, O., Toshito, T., and Uetake, T.: Study of electron identification in a few GeV region by an emulsion cloud chamber, Rev. Sci. Instrum., 74, 53,, 2003. 

Koyama, M.: Photographs of the Izu-peninsula by UAV, Iwanami Science Library 268, Iwanami Shoten, Publishers, Tokyo, 159 pp., ISBN 9784000296687, 2017 (in Japanese). 

Koyano, Y., Hayakawa, Y., and Machida, H.: The eruption of Omuroyama in the Higashi Izu monogenetic volcano field, J. Geog. (Chigaku Zasshi), 105, 475–484,, 1996 (in Japanese). 

Kreslo, I., Cozzi, M., Ereditato, A., Hess, M., Knuesel, J., Laktineh, I., Messina, M., Moser, U., Pistillo, C., and Pretzl, K.: High-speed analysis of nuclear emulsion films with the use of dry objective lenses, J. Instrum., 3, P04006,, 2008. 

Morishima, K. and Nakano, T.: Development of a new automatic nuclear emulsion scanning system, S-UTS, with continuous 3D tomographic image read-out, J. Instrum., 5, P04011,, 2010. 

Morishima, K., Kuno, M., Nishio, A., et al.: Discovery of a big void in Khufu's Pyramid by observation of cosmic-ray muons, Nature, 552, 386–390,, 2017. 

Nagahara, S. and Miyamoto, S.: Feasibility of three-dimensional density tomography using dozens of muon radiographies and filtered back projection for volcanos, Geosci. Instrum. Method. Data Syst., 7, 307–316,, 2018. 

Nakamura, T., Ariga, A., Ban, T., et al.: The OPERA film: New nuclear emulsion for large-scale, high-precision experiments, Nucl. Instrum. Methods Phys. Res. A, 556, 80–86,, 2005. 

Nishio, A., Morishima, K., Kuwabara, K., Yoshida, T., Funakubo, T., Kitagawa, N., Kuno, M., Manabe, Y., and Nakamura, M.: Nuclear emulsion with excellent long-term stability developed for cosmic-ray imaging, Nucl. Instr. Method. Phys. Res. A, 966, 163850,, 2020. 

Nishiyama, R., Miyamoto, S., and Naganawa, N.: Experimental study of source of background noise in muon radiography using emulsion film detectors, Geosci. Instrum. Method. Data Syst., 3, 29–39,, 2014. 

Nishiyama, R., Taketa, A., Miyamoto, S., and Kasahara, K.: Monte Carlo simulation for background study of geophysical inspection with cosmic-ray muons, Geophys. J. Int., 206, 1039–1050,, 2016. 

Nishiyama, R., Miyamoto, S., Okubo, S., Oshima, H., and Maekawa, T.: 3D density modeling with gravity and muon-radiographic observations in Showa–Shinzan lava dome, Usu, Japan, Pure Appl. Geophys., 174, 1061–1070,, 2017. 

Nishiyama, R., Miyamoto, S., and Nagahara, S.: Estimation of the bulk density of the Omuro scoria cone (eastern Izu, Japan) from gravity survey, Bull. Earthq. Res. Inst. Univ. Tokyo, 95, 1–7, (last access: 11 March 2022), 2020. 

Oláh, L., Tanaka, H. K. M., Ohminato, T., and Varga, D.: High-definition and low-noise muography of the Sakurajima volcano with gaseous tracking detectors, Sci. Rep., 8, 3207,, 2018. 

Rosas-Carbajal, M., Jourde, K., Marteau, J., Deroussi, S., Komorowski, J.-C., and Gibert, D.: Three-dimensional density structure of La Soufrière de Guadeloupe lava dome from simultaneous muon radiographies and gravity data, Geophys. Res. Lett., 44, 6743–6751,, 2017. 

Saito, T., Takahashi, S., and Wada, H.: 14C ages of Omuroyama Volcano, Izu Peninsula, Bull. Volcanol. Soc. Jpn. , 48, 215–219,, 2003 (in Japanese). 

Saracino, G., Amato, L., Ambrosino, F., Antonucci, G., Bonechi, L., Cimmino, L., Consiglio, L., D'Alessandro, R., De Luzio, E., Minin, G., Noli, P., Scognamiglio, L., Strolin, P., and Varriale, A.: Imaging of underground cavities with cosmic-ray muons from observations at Mt. Echia (Naples), Sci. Rep., 7, 1181,, 2017.  

Taha, A. A. and Mohamed, A. A.: Chemical, physical and geotechnical properties comparison between scoria and pumice deposits in Dhamar–Rada Volcanic Field, SW Yemen, Aust. J. Basic & Appl. Sci., 7, 116–124, (last access: 11 March 2022), 2013. 

Takahashi, S., Aoki, S., Kamada, K., Mizutani, S., Nakagawa, R., Ozaki, K., and Rokujo, H.: GRAINE project: The first balloon-borne, emulsion gamma-ray telescope experiment, Prog. Theor. Exp. Phys., Volume 2015, 043H01,, 2015. 

Tanabashi, M., Hagiwara, K., Hikasa, K., et al.: Particle Data Group, Phys. Rev. D, 98, 030001,, 2018. 

Tanaka, H. K. M., Nakano, T., Takahashi, S., Yoshida, J., Takeo, M., Oikawa, J., Ohminato, T., Aoki, Y., Koyama, E., Tsuji, H., and Niwa, K.: High resolution imaging in the inhomogeneous crust with cosmic-ray muon radiography: The density structure below the volcanic crater floor of Mt. Asama, Japan, Earth Planet. Sc. Lett., 263, 104–113,, 2007. 

Tanaka, H. K. M., Uchida, T., Tanaka, M., Shinohara, H., and Hideaki, T.: Cosmic-ray muon imaging of magma in a conduit: Degassing process of Satsuma–Iwojima Volcano, Japan, Geophys. Res. Lett., 36, L01304,, 2009. 

Tanaka, H. K. M., Taira, H., Uchida, T., Tanaka, M., Takeo, M., Ohiminato, T., and Tsuji, H.: Three-dimensional computational axial tomography scan of a volcano with cosmic ray muon radiography, J. Geophys. Res.-Sol. Ea., 115, B12332,, 2010. 

Tioukov, V., Alexandrov, A., Bozza, C., Consiglio, L., D'Ambrosio, N., De Lellis, G., De Sio, C., Giudicepietro, F., Macedonio, G., Miyamoto, S., Nishiyama, R., Orazi, M., Peluso, R., Sheshukov, A., Sirignano, C., Stellacci, S. M., Strolin, P., and Tanaka, H. K. M.: First muography of Stromboli volcano, Sci. Rep., 9, 6695,, 2019. 

Watanabe, A., Takenaka, H., Fujii, Y., and Fujiwara, H.: Seismometer azimuth measurement at K-NET station (2): Oita Prefecture, Zisin, 53, 185–192,, 2000 (in Japanese). 

Yamamoto, H.: The mode of lava outflow from cinder cones in the Ojika–Jima monogenetic volcano group, Bull. Volcanol. Soc. Jpn., 48, 1,, 2003 (in Japanese). 

Yoshimoto, M., Nakano, T., Komatani, R., and Kawahara, H.: Hyper-track selector nuclear emulsion readout system aimed at scanning an area of one thousand square meters, Prog. Theor. Exp. Phys., 2017, 103H01,, 2017. 

Short summary
In recent years, imaging of density inside volcanoes using cosmic-ray muons has been studied. In some previous studies, observations were conducted in 2 or 3 directions to obtain 3D resolution. For higher resolution, we tried to observe a volcano from 11 directions using special photographic films. The observation and analytical techniques developed in this study may be applied to other volcanoes and large objects. Reconstructing 3D images of the volcano is in progress using the obtained data.