Three-dimensional density tomography determined from multi-directional muography of the Omuroyama scoria cone, Higashi–Izu monogenetic volcano field, Japan

Scoria cone is one of the basic forms of volcanoes, and clarifying its detailed internal structure and magma movement during its formation is important for understanding the nature of eruptions and also for volcanic disaster mitigation. We conducted a multi-directional muographic survey of the Omuroyama scoria cone, Japan, in order to determine the three-dimensional density structure of the scoria cone. We used a nuclear emulsion detector optimized for multi-directional muography that was installed at 11 sites around the volcano. Muon tracks recorded on emulsion films were read with a high-speed automated readout system. We obtained the three-dimensional density structure by applying a linear inversion to the muographic images, and evaluated the uncertainties. High-density zones were detected in scoria cone. Based on these observations and detailed topographic and geological constraints, we infer that high-density zones are the central, highly welded vent of the scoria cone, three-directional radial dikes extending from the central vent, respectively. We also infer that an E-W-trending dike fed a small lava flow, and that a SSE-NNW-trending dike caused a small flank eruption and produced a crater on the flank of the cone. Our results visualize the three-dimensional internal structure of volcanoes with novel resolution, and thereby demonstrate that the formation process of volcanoes can be discussed in detail. Since it has been difficult to directly detect such clear visualization of the internal structures and eruptive processes in the past, multi-directional muography has resulted in providing one hope for understanding volcanic activity.


Introduction
Scoria cones are a type of monogenetic volcano that occurs worldwide. Previous studies of scoria cones have revealed a variety of internal structures and developmental processes, and the factors responsible for these features are partly understood (McGetchin et al. 1974;Yamamoto 2003;Kereszturi and Nemeth 2013;Smith and Nemeth 2017). However, these studies are mostly based on fragmented outcrops and remnants of cones, where much of the edifice has already been eroded. There have been few case studies of cross-sections of scoria cones from multiple directions that can examine the relationships with the dikes and lava flows associated with the scoria cone. One of the most significant examples is a study that clarified the three-dimensional relationships between the central conduit of the Puy de Lemptegy scoria cone and the dikes and lava flows derived from it, including the direction of magma movement, were determined from a well-exposed and well-preserved quarry site (Petronis et al. 2013;Delcamp et al. 2014). However, this is a rare case due to the small size of the scoria cone and the fortunate fact that the entire quarry site was preserved without planting. The detailed, three-dimensional, inner structures of large scoria cones (i.e., edifices with talus slopes classified as developmental stages 3-4 in McGetchin et al. (1974)) have not been determined.
A few studies have conducted indirect geophysical exploration of interior of scoria cones. Hashimoto et al. (2009) obtained the resistivity and magnetic structures of Komezuka scoria cone on Aso Volcano by electromagnetic surveys. This study identified a high-density area in the center of the scoria cone, which was presumed to be welded, and showed that the central conduit was not connected to the surrounding lava flow, at least in the shallow part of the cone. Barde-Cabusson et al. (2013) conducted an electrical survey of three scoria cones in Spain and obtained results that identified a central conduit and surrounding welded rocks. Bolós et al. (2020) summarized the results of a selfpotential, ground-temperature, and a 3D resistivity model geophysical survey on the Parícutin scoria cone, and they understood the surface structure and detected zones of convective hydrothermal fluid circulation in the subsurface. However, in these cases, the inner structure of the cones was only roughly estimated by comparison with the structures observed in geological surveys of other volcanoes.
Muography has rapidly developed since the twenty-first century as an exploration method that can reveal the internal structures of large objects (Tanaka et al. 2007a). This technique uses the absorption of cosmic ray muons in matter to measure the density distribution like an X-ray image. The attenuation of the muon flux depends on the amount of material along the path of the muon (i.e., the density length). The energy spectrum of the incident cosmic ray muon (Honda et al. 2004) and its interaction with matter (Groom et al. 2001) allows the density length to be obtained from measurements of the penetrating muon flux.
Unidirectional muography has been conducted at numerous volcanoes (Mt. Asama; Tanaka et al. 2007b; Showa-Shinzan lava dome; Tanaka et al. 2007c;Satsuma-Iwojima;Tanaka et al. 2009Tanaka et al. , 2014Sakurajima;Oláh et al. 2018;Puy de Dôme;Carloganu et al. 2013;Ambrosino et al. 2015;Stromboli;Tioukov et al. 2019). These observations only obtained two-dimensional projection images like radiographs. To solve this problem, Tanaka et al. (2010) undertook muography from two directions. Other studies also used a combination of muography and gravity observations, which are also sensitive to the density (Nishiyama et al. 2014a(Nishiyama et al. , 2017Rosas-Carbajal et al. 2017). However, these approaches were poor three-dimensional resolution. For example, the resolution of Rosas-Carbajal et al. 2017 is approximately 500 m. The simulation of Nagahara and Miyamoto (2018) showed that in order to solve this problem and obtain an appropriate three-dimensional density distribution, it is necessary to conduct muographic observations in multiple directions, as the case for X-ray computed tomography. For example, they showed that observing from 64 directions would provide a spatial resolution of 10 to 20 meters.
Such multi-directional observations are difficult due to detector limitations, but have become feasible in recent years with the rapid development of muon detectors. In order to perform multi-directional muographic observations of volcanoes, it is necessary to install the detectors in places where infrastructure is not available. Nuclear emulsion muon detectors are compact, lightweight, and do not require a power supply, which makes these suitable for multi-directional observations. With the development of related technologies in recent years (Hamada et al. 2012;Yoshimoto et al. 2017), emulsion detectors have also become increasingly larger in scale (for example, Khufu's Pyramid; Morishima et al. 2017).
Scoria cones are suitable targets for multi-directional muography. Given that the interiors of scoria cones are expected to have a complex structure due to diking and welding (Kereszturi and Németh 2016). It is challenging to explore their structures with conventional muographic observations. However, multi-directional muography might allow scoria cones with complex structures to be robustly investigated.
We selected the Omuroyama scoria cone (34° 54′ 11′′ N, 139° 05′ 40′′ E; 580 m above sea level.), which is the largest monogenetic scoria cone in the Higashi-Izu monogenetic volcano field located in Shizuoka Prefecture, Japan (Koyama and Umino 1991), for a multi-directional muographic study (Fig. 1). The diameter of the foot of the edifice is about 1 km, and the height from the bottom is 300 m. There is a summit crater with a diameter of ca. 300 m and a depth of ca. 50 m. The summit crater is surrounded by straight slopes that strongly suggest its talus origin. Saito et al. (2003) estimated the eruption age of the cone to be ~ 4 ka, based on radiocarbon dating.
The eruptive history of the Omuroyama scoria cone has generally been constrained by the previous field geological survey (Koyano et al. 1996). The process of the eruption of the volcano can be divided into five stages (A to E), based on tephra stratigraphy. In addition, associated lava flows were divided into lava flows I to IV. The eruptive history is as follows.
Eruption of the Omuroyama scoria cone began with a subplinian eruption (stage A) and a short period of deposition of volcanic ash (stage B). Strombolian eruptions then continued for a protracted period, which formed a large scoria cone and lava flows I-III from three locations at the foot of the cone (stages C-D). Lava flows II and III formed the Izu Kogen lava plateau and Jogasaki lava coast, respectively, and finally the Iwamuroyama and Moriyama plug domes formed on the source craters of lava flows II and III, respectively. At the end of stage D, a lava lake was formed in the summit crater, which then drained. Part of the drained lava leaked into the western foot of the cone, which produced a small-scale lava flow (lava flow IV). Finally, the volcanic activity ended with a small explosion on the southern flank of the cone, which formed a small crater (stage E).
Because no strata showing significant dormant periods of eruptive activity have been identified in the tephra succession of the Omuroyama scoria cone, Koyano et al. (1996) also concluded that the Omuroyama scoria cone should be a monogenetic volcano and its eruption time was estimated to be a few years, 20 to 30 years at the most.
Although the Omuroyama scoria cone has an apparently well-proportioned cone-shaped edifice, a small ridge  Koyama) on its western foot, a small crater on the southern flank, and the complex topography of the summit crater imply a complex formation history and inner structures (Fig. 1). In typical basaltic rocks, the density of the non-welded scoria deposits is expected to be 0.5-1.0 × 10 3 kg/m 3 , whereas the density of the lava flows and welded scoria are expected to be > 2.0 × 10 3 kg/m 3 . With this amount of density contrast, we can expect to obtain a resolution of less than 100 m by muography from a dozen or more directions (Nagahara and Miyamoto, 2018). In addition, no mountains, which would interfere with the muographic observations, exist behind the Omuroyama scoria cone in any direction within a radius of 4 km. Therefore, the Omuroyama scoria cone is an ideal target for multi-directional muography.
We undertook multi-directional muographic observations of the Omuroyama scoria cone and used the data to obtain a three-dimensional density reconstruction using linear inversion techniques. In addition, we used our reconstruction, the latest topographic data, and field geological surveys to constrain the eruption history of the Omuroyama scoria cone since its formation.

Overview of multi-directional muography
The multi-directional muography and analysis of the results can be divided into several steps. Figure 2 shows the steps taken when using emulsion detector and related previous studies. A feasibility study on Omuroyama is described in Nagahara and Miyamoto (2018). The procedure to obtain the two-dimensional muographic images is described in Miyamoto et al. (2022). Three-dimensional density reconstruction and its interpretation is described in a later section of this paper.

Muographic observations
The muographic observations were performed in 2018 and 2019. The detail of muon detector, the installation of detectors, and analysis are described in Miyamoto et al. (2022). Figure 1 shows the installed stations. To calculate the path length, we used a 1-m resolution digital elevation model (DEM) using an airborne LiDAR survey obtained in 2009 by the Numazu Office of River and National Highway, Ministry of Land, Infrastructure, Transport and Tourism.
The detector was newly constructed based on the structure of an emulsion cloud chamber (ECC; Nishiyama et al. 2014b). Because the background noise that affects muographic studies is due to protons, electrons, and muons with energies of < 1 GeV (Nishiyama et al. 2016;Niess, 2022), 1-mm-thick lead plates were stacked between the nuclear emulsions. In this study, 20 emulsions with sizes of 10 cm × 12.5 cm and 19 lead plates of the same size were alternately stacked .
After the observations, the nuclear emulsion films were developed and the recorded tracks were read using the hypertrack selector (HTS; Yoshimoto et al. 2017). Reconstruction of the tracks was performed with NETSCAN 2.0 software (Hamada et al. 2012). We then selected the muon tracks and estimated the muon detection efficiency .
The ratio of the observed/expected muon counts in the free sky region was calculated for each 100-mrad angular bin, and the mean value was used as the calibration value and the deviation as the calibration error (Table 1). Only ECC_C (station NNE) had a larger calibration error than the other detectors. We subsequently analyzed the data from 10 of the 11 sites, excluding this station.
We calculated the density length using the following procedure.
1. We calculated a numerical table of the relationship between the mean density and muon counts for each angular bin using the energy spectrum of the incident cosmic ray muons from Honda et al. (2004), their interaction with matter from Groom et al. (2001), and topographic data for the Omuroyama scoria cone. 2. We multiplied the expected muon counts by the calibration value of the ECC. 3. We determined the density length that matched the observed muon counts. The feasibility study is shown in detail in Nagahara and Miyamoto (2018). The obtaining of two-dimensional muographic images is shown in detail in Miyamoto et al. (2022), track scanning is shown in Yoshimoto et al. (2017), and track reconstruction is shown in Hamada et al. (2012), respectively An example of the two-dimensional muographic results is shown in Fig. 3. The overall density of the Omuroyama scoria cone was 1.39 ± 0.01 × 10 3 kg/m 3 . This is in good agreement with the value obtained by Nishiyama et al. (2020) (1.39 ± 0.07 × 10 3 kg/m 3 ) from gravity observations.

Inversion method
The three-dimensional density reconstruction was undertaken using the observed muon counts based on the methodology presented in this section.

Three-dimensional density inversion method
We used the linear inversion method of Nishiyama et al. (2014a). The density length d can be written as follows: where ρ is the density, τ is the muon path, and ξ is the coordinate along the trajectory τ. The region of reconstruction is subdivided into J voxels of density j ( j = 1, 2, … , J ). The density length of path i ( i = 1, 2, … , I ) can be written as follows: where d i is density length of path i and A ij is the length of the ith path across the jth voxel. Equation 2 can also be written as follows: where d is the data vector, is the observation matrix, and is the voxel density. Nishiyama et al. (2014a) used a Bayesian approach to solve this inverse problem (Tarantola and (1)  Nercessian, 1984). When the observed data d is obtained, the solution of Eq. 3 is as follows: where 0 is the initially guessed density, d is the data covariance matrix, is a covariance matrix for the initial guess, and ′ is posterior covariance matrix. Nishiyama et al. (2014a) imposed a smoothing constraint through an exponential covariance on , based on Tarantola and Nercessian (1984): where is the magnitude of the uncertainty, l is the correlation length, and D j 1 , j 2 is the distance between the j 1 th and j 2 th voxels. Procedures for determining the optional parameters, l and , are described in later sections.
In this method, the muon arrival direction is represented by the center of the angular bin because density length cannot be defined for a finite range of solid angles. This is not a good approximation when the angular resolution of the muographic data is poor. To solve this problem, we defined the observation matrix in terms of the number of observed muons per angular bin. The ith observed muon count is the sum of the muon numbers in the smaller angular bins i . Therefore, we first prepared an observation matrix with the observed muon counts for i , and then summed over to calculate the observation matrix for the observed muon counts of the solid angle Ω i (Fig. 4).
The relationship between the cosmic ray muon flux attenuation and density length is nonlinear (Tanaka et al. 2007b). Therefore, the density distribution within the scoria cone is linearized here for the case of 0 . When the relationship between the muon number N and density length d is N Obs = f (d) , the density length when the density distribution is 0 is d 0 , the muon count is N 0 Obs , and the observed muon count in the angular bin i is as follows: Using this and Eq. 2, the observation matrix for the number of observed muons A N can be expressed as follows: Summing this over , we obtain the observation matrix for the solid angle Ω i : In this case, the data covariance matrix for the muon counts is as follows: where 2 N i is the error of N i for the i-th solid angle area. From the above, the formula for the three density reconstruction is: The calculations for this analysis can be done on a laptop, but multiple calculations are required to determine optional parameters. The method of determining the parameters in this study is described in a later section.

Topographic and geological survey
To interpret the muographic results, we undertook new detailed topographical and geological surveys of the Omuroyama scoria cone and surrounding region. For the topographic analysis, we created a high-resolution DEM (10-cm mesh), based on the point cloud data acquired by Shizuoka Prefecture using an airborne LiDAR survey in 2019 (https:// www. geosp atial. jp/ ckan/ datas et/ shizu oka-2019-point cloud). A field geological survey was also conducted from January to June 2021 to interpret the geomorphic features. Their results are shown in the later section.

Three-dimensional density reconstruction
The inversion method was applied to the muographic results from 10 directions and 3D density reconstruction was performed. In this subsection, we describe the determination of the reconstruction parameters ( l and ) using the observed data, the resulting estimation of uncertainties, and the results of the resolution estimation.

Parameter determination
For the three-dimensional density reconstruction, we selected the following parameters: voxel size of 20 × 20 × 20 m, 0 = 1.39 × 10 3 kg/m 3 , = 0.3 × 10 3 kg/m 3 , and l = 70 m. For 0 , we used the observed value of the overall density. The above values were selected by searching for parameters for l in increments of 10 m from 10 to 200 m, and for in increments of 0.1 × 10 3 kg/m 3 from 0.1 to 1.0 × 10 3 kg/m 3 . The three-dimensional density reconstruction results for each parameter are shown in Fig. 5.
The correlation length was determined by requiring that the amount of information in the three-dimensional density reconstruction results and observed data be equal. Therefore, the amount of information in the three-dimensional density reconstruction results was considered to be the number of voxels when the length of one side of the voxel was the correlation length. Given that the correlation length at which the amount of information in the three-dimensional density reconstruction is equal that in the observed data is 63.6 m, we selected 70 m.
We constrained the value of in two ways. The first constraint is that the results of the three-dimensional density reconstruction should have realistic values. In this case, the constraint was that > 90% of the voxels should have densities of 0.5-3.0 × 10 3 kg/m 3 . As a result, = 0.1-0.4 × 10 3 kg/ m 3 was applicable at l = 70 m. The second constraint is that the structure is reconstructed by simulations assuming a 100-m diameter, 2.5 × 10 3 kg/m 3 density cylinder structure on a mountain body of 1.3 × 10 3 kg/m 3 density. The simulation results for the 450-m-elevation section are shown in  6. As a result, the cylindrical structure was not reconstructed for = 0.1 × 10 3 kg/m 3 . From these considerations, we conclude that = 0.2-0.4 × 10 3 kg/m 3 . In this range, the density contrast changed, but there was no significant difference in the density structure. Therefore, we selected the central value of 0.3 × 10 3 kg/m 3 .

Uncertainties
The three-dimensional density reconstruction results include a variety of uncertainties. In this paper, the following three uncertainties were considered.
The first uncertainty is the covariance matrix. The propagation of the accidental error on the observed data is represented by the diagonal component of the covariance matrix, which is represented by Eq. 11. The resulting uncertainty was 0.2-0.3 × 10 3 kg/m 3 throughout the scoria cone.
The second uncertainty is 0 . In this study, we used 0 as the observed value of the overall density. There is an uncertainty on 0 that can be assigned. In this case, the variation of the overall density determined for each observation point was taken into account. The difference between the results of the three-dimensional density reconstruction for 0 = 1.35 and 1.44 × 10 3 kg/m 3 , and those for 0 = 1.39 × 10 3 kg/m 3 was adopted as the second uncertainty. The uncertainty was less than 0.1 × 10 3 kg/m 3 and smaller than that of the covariance matrix.
The third uncertainty is the stability of the solution depending on location of the observation points. In a threedimensional density reconstruction obtained by an inversion method, some artificial noise can characterize locations with little observational information. In order to evaluate the effect of artificial noise related to the station configuration, we performed three-dimensional density reconstructions using 9 observation points by removing 1 of the 10 observation points, and examined the changes in the density reconstructions. The results are shown in Fig. 7. It was found that high-density regions tend to form in front of the removed observation point and low-density regions form on both sides of the observation point. In this study, the uncertainty of each voxel was defined by the value of the largest density change of these 10 analysis results.
The three uncertainties were assumed to be independent. The square root of the sum of the squares of the above three uncertainties was taken as the overall uncertainty of the three-dimensional density reconstructions, and the results and topography from the 10 cm DEM are shown in Fig. 8. The high-density structures identified by the muography Fig. 6 A simulation in the case where l = 70-100 m and = 0.1-0.5 × 10 3 kg/ m 3 , which is at an elevation of 450 m. The structure was assumed to be cylindrical with a diameter of 100 m and density of 2.5 × 10 3 kg/m 3 below z = 500 m, with a density of 1.3 × 10 3 kg/m 3 . The blue dots represent the locations of the observation points, and the black lines represent the contours of the scoria cone in the cross-section observations can be classified into five parts: (1) the central part of the edifice just below the summit crater; (2A-C) the westward, south-southeastward, and north-northeastward extensions from the central part of the edifice; (3) the southwestern edge of the Iwamuroyama lava dome. Figure 9 shows a view of the three-dimensional density image of the Omuroyama scoria cone obtained by multi-directional muography and high density parts. Interpretation of these structures is given in later sections.

Evaluation of the resolution
In order to assess the resolution of the reconstruction results, we simulated the following two cases of simple structures: (a) a scoria cone density of 1.3 × 10 3 kg/m 3 and density of a 50-m-thick dike is 2.5 × 10 3 kg/m 3 ; (b) a scoria cone density of 1.3 × 10 3 kg/m 3 and density of a cylindrical structure with a diameter of 100 m is 2.5 × 10 3 kg/m 3 ; (c) a scoria cone density of 1.3 × 10 3 kg/m 3 and density of two cylindrical structures 250 m distant from each other with a diameter of 100 m is 2.5 × 10 3 kg/m 3 . The results are shown in Fig. 10. Two main features were identified. Firstly, the thickness of the dike and radius of the cylinder tend to be blurred over a thickness that is twice the size of these bodies. Secondly, low-density artificial noise occurs around the high-density areas. It was also identified that the two columns, 250 m distant from each other, were clearly distinguished. This should be taken into consideration when interpreting the results of the three-dimensional density reconstructions.

Topographic and geological results
High-resolution topographical map and its interpretation of each topographic feature is shown in Figs. 1b and 11. The ridge at the western foot of the Omuroyama scoria cone (Fig. 1b) was confirmed to be lava flow IV, which emanated from the base of the scoria cone at the end of the eruption (Koyano et al. 1996). We collected samples of this lava flow and measured their density, which was 2.22 × 10 3 kg/m 3 .
The small craters located on the south flank of the scoria cone have a double crater rim (Fig. 1b), indicating that the crater position shifted slightly during the eruption. The shallow trench-like depression extending from the summit crater rim to the south craters is thought to have been caused by slope collapse of surface deposits during the flank eruption.
The topography of the summit crater is rather complex, reflecting the final stage of the eruption (Fig. 11c, d). There is a flat-bottomed depression (i.e., the crater floor). In the outcrop on the northeastern inner wall of this depression, an undeformed spindle-shaped volcanic projectile can be observed in the dense lava flow, which is interpreted to have fallen into the lava lake that previously existed in the crater (Koyano et al. 1996). This means that a lava lake with a level of at least 535 m above sea level existed during the eruption and was then drained at the end of the eruption. Note that the high-density area that may correspond to this lava lake Fig. 8 The three-dimensional density reconstruction results, topography, and estimated overall uncertainties of each elevation section. The blue dots represent the locations of the observation points. (i) Three-dimensional density reconstruction results shown with 50% transmission in the images. (ii) The overall uncertainties fragment is not visible (Fig. 8), suggesting that this fragment is a minor structure and that the majority of the lava lake has drained back.
There are several ridges (northern, mid-western, and northeastern) that appear to reflect strongly welded zones near the surface. A small depression (Kitsune-no-oyado), which occurs near the northern edge of the crater rim, is thought to be a very small crater produced by a small explosive event. There is also a small cinder cone in the northwestern bottom of the summit crater.
Northeast of the summit crater, a pit that was 1 m in diameter and 4 m deep was excavated in 2016 (Fig. 12).

Fig. 9
View of a three-dimensional density image of the Omuroyama scoria cone from a the sky above the crater, b the eastern sky, and c) the southwestern sky. The volumes of high density are shown with high opacity. The elevation is scaled to be double that of the hori-zontal scale. The orange dots represent the locations of the observation points. The green numbers mark the high-density structures. This figure is also viewable at https:// www. eri.u-tokyo. ac. jp/ CHEER/ data/ omuro 3ds/ This site is significant because of the lack of large outcrops on the surface of the Omuroyama scoria cone. Most of the outcrop is covered by red and strongly oxidized scoria, with some gray, strongly welded rocks. Samples were collected at depth intervals of 1 m from the pit, which have densities of 1.32-2.36 × 10 3 kg/m 3 . In general, the density increased with depth. These value are consistent with the value of the three-dimensional density reconstruction.

Structure and formation of the Omuroyama scoria cone
In this section, we interpret the results of the three-dimensional density reconstructions (Fig. 9) and discuss the formation of the Omuroyama scoria cone, based on the results of the detailed topographic and geological surveys. The central part of the edifice just below the summit crater has a vertical and cylindrical shape and is inferred to be the main vent, which is highly welded. The 2A extension propagates horizontally from the base of the main vent towards the source of lava flow IV. As such, 2A can be interpreted as being a dike that fed lava flow IV. This implies that draining of the lava lake in the summit crater was caused by the dike intrusion that formed lava flow IV.  The 2B extension also propagates horizontally from the base of the main vent towards the SSW. This feature is also considered to be a dike, but there is no geological evidence for this. However, the small craters on the southern flank might have been supplied by this dike.
The 2C extension is shorter and less clear than the other high-density features. Moreover, it is possible this represents artificial noise, because this is a location where there is no observation point. However, we interpret 2C to represent an incipient dike because of the triaxial elongation of 2A-C, because such a triaxial intrusive structure would be expected to result from expansion of the volcanic edifice caused by magma input during the eruption.
Similar three-directional dike systems have been observed in other volcanoes worldwide, including at Shiprock (Navajo Volcanic Field, northwestern New Mexico: Townsend et al. 2015), Hawaii (Wyss 1980), and the Canary Islands (Carracedo and Troll 2013). The Higashi-Izu monogenetic volcano field, in which the Omuroyama scoria cone occurs, is located in a NW-SE compressional stress field caused by the collision between the Izu and Honshu arcs, and has vent alignments of the same azimuth (Koyama and Umino 1991). This situation appears to be inconsistent with the radial orientation of the dikes associated with the eruption of the Omuroyama scoria cone. While volcanic edifices without significant relief tend to develop subparallel dike swarms consistent with the regional tectonic trend, edifices with subcircular topographic relief are characterized by radial dikes under limited regional stress control (Acocella and Neri 2008). Since the Omuroyama scoria cone is an isolated mountain, its interior and surroundings are less susceptible to the influence of the regional stress field and thus was probably affected by the local stress field caused by magma ascent, producing the radial dikes.
Although the high-density feature (3) is associated with a large error (> 1.0 × 10 3 kg/m 3 ), this structure is consistent with the Iwamuroyama plug dome that formed at the end of the lava eruptions due to the increased lava viscosity.
Based on the above, a tentative history of the formation of the Omuroyama scoria cone is proposed (stages C-E outlined in introduction). Stages A and B that occurred prior to construction of the scoria cone were not considered. Most of the scoria cone was formed during stage C when lava flow I was also effused in the early part of this stage (Fig. 13a). In stage D, lava flows II and III were effused from the northeastern and southern foot of the volcano, respectively (Fig. 13b). Both lava flows are distant from the center of the Omuroyama scoria cone, unlike lava flow I. Therefore, it is reasonable to assume that the lava flows were fed by a twodirectional radial dike intrusion from the central vent. This dike intrusion would have been sourced from the deeper part of the Omuroyama central vent, but the latter is not visible in our muographic observations. Two lava domes (i.e., plug domes), Iwamuroyama and Moriyama, were formed at the outlets of lava flows II and III, respectively, due to the increase in lava viscosity when the lava flows stopped. As a result, dike penetration ceased, and the magma ascended to form a lava lake in the summit crater (Fig. 13c).
In stage E, a further increase in pressure caused radial dikes to intrude from the central conduit in three directions, which were to the NNE, W, and SSE (Fig. 13d). The W-propagating dike formed lava flow IV at the western foot of the cone. The SSE-propagating dike formed a small crater on the southern flank of the come. The decrease in pressure caused by the formation of these dikes led to the lava lake in the summit crater becoming drained (Fig. 13e).
Even after the lava lake drained, the area around the central conduit remained hot and welding proceeded to form a high-density area in the central part of the edifice. The remaining magma and volcanic gases formed small pyroclastic deposits and craters in the western part of the summit.

Possibility of applying multi-directional muography to disaster mitigation
Let us consider possible applications of multi-directional muography to volcano disaster mitigation. One main result of this study was the visualization of the relationship between the formation of dikes within the volcanic edifice and eruptive processes of the scoria cone, particularly the lava flows at the foot of the cone and formation of flank craters.
Lava flows from the base of volcanos, not only scoria cone, can be a major risk to residential areas, and such flows can lead to collapse of the volcanic edifice. In a scenario focused on the summit eruption, the sudden occurrence of flank eruptions can also be significant. Therefore, it is important for volcanic disaster mitigation to forecast the phenomena that may occur on the flank and surrounding the volcano during an eruption. However, it has been difficult to detect magma movement occurring inside such volcanos. To what extent can multi-directional muography contribute to solving these problems?
First, we discuss the spatial resolution of muography. The Omuroyama scoria cone is a rare example where muon detectors can be installed around the volcanic edifice from all directions. In most cases, the topography, vegetation, and access limit the observation directions. Therefore, in the case of real volcanoes, it is difficult to select the location of observation points that can provide the required resolution. When the number of observation directions is small, the resolution is low (Nagahara and Miyamoto 2018). Observation direction bias causes density bias as seen in Fig. 7. To overcome such cases, it is required to develop an analytical method to obtain three-dimensional density reconstruction that reduces the bias even when the location of observation points is not ideal.
Next, it is necessary to mention the temporal resolution that current muographic techniques have reached: Tanaka et al. 2014 succeeded in detecting the ascending and descending of the magma head of Satsuma Iwo Jima volcano in Japan as time-lapse movies taken every 3 days. Also, Jourde et al. (2016) successfully detected fluid movement due to hydrothermal system activity at La Soufriere volcano, Guadeloupe, France, with a temporal resolution of 1 week. Current techniques are still inadequate because the time scale for lava to flow outward from a scoria hill due to dyke intrusion is often shorter than 1 day. Improving the temporal resolution of muography also requires increasing the number of detected muons that penetrate the mountain body.
On the other hand, the flux of cosmic ray muons is limited, the area of the sensitive detection surface must be increased. Also, in order to increase the effective solid angle, the detector must be placed as close to the volcano as possible. However, it is usually difficult to approach an active volcano. Therefore, it is necessary to develop techniques for remote installation of muon detectors. In addition, multi-directional observations require more muon detectors than in the previous studies described above. Therefore, making 3D density image with higher temporal resolution than these previous studies by multi-directional muography still has a lot of technical challenges to be solved.

Possibility of applying multi-directional muography to other targets
Let us consider about the application of multi-directional muography to other fields. As described in section evaluation of the resolution, in Omuroyama, we were able to detect structures larger more than 100 m for a density difference of about 2 times. Therefore, objects that are similar in size and density difference to the vent of Omuroyama are expected to be detectable. Such examples include volcanoes as well Fig. 13 The five major stages of development of the Omuroyama scoria cone. a Stage C when most of the scoria cone was formed and lava flow I was effused. b Stage D when lava flows II and III were effused. c Lava effusion ceased and the magma ascended to form a lava lake in the summit crater. d Stage E when a further pressure build-up caused radial dikes to intrude outward from the central conduit in three directions. e The W-propagating dikes that formed lava flow IV. The SSE-propagating dike caused an eruption and formation of the small craters on the southern flank of the scoria cone. The pressure decrease due to the formation of the dikes caused the lava lake in the summit crater to drain as archeological sites (Khufu's Pyramid; Morishima et al. 2017) and nuclear reactors (Fukushima Daiichi; Fujii et al. 2022).
To detect smaller structures, it is required to increase the number of observation directions since the detectable structure size strongly depends on it (Nagahara and Miyamoto, 2018). For example, to halve the diameter of the detectable structure, it is required to increase the number of observation directions by four times.
On the other hand, when density differences are small, the uncertainty of the image should be reduced. In this case, not only should the number of observation directions be increased, but the error in the muographic image at each observation point should be reduced. This requires increasing the number of muons observed at each station. The accidental error in muography at each station depends on the statistics of the muon. For example, to halve this error, it is required to increase the muon statistics by four times.

Conclusions
Muographic observations were obtained from 11 directions of the Omuroyama scoria cone, and three-dimensional density reconstructions were made of the data collected from 10 of the 11 directions that were successfully observed. A high-density cylindrical region exists just below the summit crater in the central part of the scoria cone, from which two prominent high-density zones extend to the W and SSE, and a small high-density zone also extends to the NNE. Based on these results and topographical and geological surveys, we constrained the formation of the Omuroyama scoria cone.
One main result of this study was the visualization of the relationship between the formation of dikes within the volcanic edifice and eruptive processes of the scoria cone, particularly the lava flows at the foot of the cone and formation of flank craters. It is important for volcanic disaster mitigation to forecast the phenomena that may occur on the flank and surrounding the volcano during an eruption. Since it has been difficult to directly detect such features in the past, multi-directional muography has resulted in providing one hope.
As mentioned in the previous section, current muography has a time resolution of only a few days, making it difficult to detect signs of a rapid change of eruption in advance even with multi-directional muography. It is also difficult in real active volcanoes to install a large number of detectors surrounding a very close location of the volcano in order to obtain high resolution. However, if more muon fluxes can be detected in a shorter time and detectors can be remotely installed at active volcanoes with improved detectors in the future, it will be possible to confirm the stability of scoria cones once formed, for example, or to detect signs of collapse in the interior of a mountain with sufficient lead time.
In addition, this study demonstrates that detailed threedimensional density reconstructions can be obtained from muographic observations in multiple directions. This approach can be widely applied in other fields, such as the investigation of archeological sites and other structures.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.