Muon Tomography sites for Colombian volcanoes

By using a very detailed simulation scheme, we have calculated the cosmic ray background flux at 13 active Colombian volcanoes and developed a methodology to identify the most convenient places for a muon telescope to study their inner structure. Our simulation scheme considers three critical factors with different spatial and time scales: the geomagnetic effects, the development of extensive air showers in the atmosphere, and the detector response at ground level. The muon energy dissipation along the path crossing the geological structure is modeled considering the losses due to ionization, and also contributions from radiative Bremstrahlung, nuclear interactions, and pair production. By considering each particular volcano topography and assuming reasonable statistics for different instrument acceptances, we obtained the muon flux crossing each structure and estimated the exposure time for our hybrid muon telescope at several points around each geological edifice. After a detailed study from the topography, we have identified the best volcano to be studied, spotted the finest points to place a muon telescope and estimated its time exposures for a significant statistics of muon flux. We have devised a mix of technical and logistic criteria --the ``rule of thumb'' criteria-- and found that only Cerro Machin, located at the Cordillera Central (4$^{\circ}$29'N 75$^{\circ}$22'W), can be feasibly studied today through muography. Cerro Negro and Chiles could be good candidates shortly


Introduction
Muography is a non-invasive astroparticle technique, introduced several decades ago for imaging anthropic and geologic structures, benefiting from the high penetration atmospheric muons produced in cosmic rays showers. This method offers a spatial resolution in the order of tens of meters, up to a kilometer of penetration, and measures the out-coming muon flux for different directions by means of a hodoscope. The flux variance between trajectories allows us to extract information about the inner density distribution of the scanned object.
In the beginning, the study of cosmic rays -and particularly their detection after crossing geological structures-were motivated by the need to understand the background noise in particle detectors inside anthropic/natural structures [1,2]. Luis Álvarez [3] was the first to apply muon radiography, with no results, to the pyramid of Cheops; but this initiative led to the successful ScanPyramids project, which recently discovered cavities in the Pyramid of Khufu (Cheops) [4].
The interest in muon radiography for Earth sciences studies arose after the discovery of the significant penetrating power of some high energy secondary particles produced by the interaction of cosmic rays with the atmosphere (hereafter, primaries). These particles can cross hundreds of meters of rock with an attenuation related to the amount of matter traversed along its trajectory [5]. This technique uses the same basic principles as a standard medical radiography: it measures, with a sensitive device, the attenuation of cosmic muons when crossing geological structures. Although there are limitations with muography in detecting deep structures beneath the volcano (such as a magma chamber), it is particularly useful when it is applied to shallow volcano phenomena such as the conduit dynamics. Thus, volcano muography constitutes an unique method to obtain direct information on the density distribution inside geological objects with a better spatial resolution than other geophysical techniques (see [6,7,8,9,10,11,12,13,14,15,16,17,18], and references therein).
Colombia, located in the Pacific Belt, has more than a dozen active volcanoes (see figure 1) clustered in three main groups along the Cordillera Central, the highest of the three branches of the Colombian Andes. Most of these volcanoes represent a significant risk to the nearby population in towns and/or cities [19,20,21] and have caused major disasters. The most recent, the Armero tragedy occurred in November 13, 1985, when pyroclastics of the Nevado del Ruiz fused about 10% of the mountain glacier, sending lahars with the terrible, devastating result of 20,000 casualties [22]. Therefore, to determine and to model the inner volcano structure is crucial in evaluating its potential risks. Muon tomography is a powerful technique, which measures the cosmic muon flux attenuation by rock volumes of different densities, allowing the projection of images of volcanic conduits at the top of the volcanic edifice. It constitutes an attractive way to infer density distributions inside different geological structures, which is critical in the study of possible eruption dynamics associated with specific eruptive styles. Nowadays, astroparticle research groups in Colombia have started to explore this technique to estimate the density distribution within geological edifices by recording the variation of the atmospheric muon flux crossing the structure [23,24,25,26,27,28,29].
The objective of the present work is twofold: first, to identify possible muon telescope observation sites to study Colombian active volcanoes by estimating the muon flux emerging from the geological edifices and second, to determine the time exposure for our hybrid muon telescope at those selected sites. Section 2 describes four of the main active Colombian volcanoes. In Section 3, we discuss the detailed simulation scheme used to estimate the incoming atmospheric muon flux and, we also investigate the energy loss of muons crossing the volcano edifice. For completeness, our hybrid muon telescope is briefly described in section 4. In 5, we analyze muon propagation across the Machín volcano and calculate the exposure times for four observation points. We also estimate the outgoing muon flux from the geological structures that could be detected by our telescope. Besides, this Section discusses a ray-tracing analysis for muon propagation at three other active volcanoes: Chiles, Cerro Negro, and Galeras. From the previous results we devise in Section 6 "rule of thumb" criteria that should be fulfilled by the potential sites (listed below) and apply them to 13 active Colombian volcanoes. Finally, in Section 7 conclusions are presented along with possible future works.

Volcanoes in Colombia
In this work, we have considered 13 Colombian active volcanoes: Azufral, Cerro Negro, Chiles, Cumbal, Doña Juana, Galeras, Machín, Nevado del Huila, Nevado del Ruíz, Nevado Santa Isabel, Nevado del Tolima, Puracé, and Sotará. Figure 1 displays an artistic representation of their geographic distribution. Because of their social significance and eruptive history, we shall briefly describe some of the characteristics of four of them: Galeras, Nevado del Ruiz, Cerro Machín, and the Cerro Negro-Chiles complex.

Cerro Machín
The Machín volcano is often overlooked as a minor edifice in the Cerro Bravo-Cerro Machín volcanic belt but, considering its high explosive potential, dacitic composition and magnitude of past eruptions, it must be considered one of the most dangerous active volcanoes in Colombia [30]. It is also located at the Cordillera Central (4 • 29'N 75 • 22'W), between Cajamarca and Machín [31].
Over the last 5,000 years, Machín has had six eruptions -the last one occurred about 850 years agogenerating pyroclastic flows, depositing several tens of centimeters of ashlayers, throwing eruptive columns (several tens of kilometers) and flows of volcanic mud. In recent times, some of the manifestations of Machín's volcanic activity are the presence of fumaroles, permanent micro-seismicity, thermal waters flowing in the vicinity of the crater, geoforms of the well-preserved volcanic building and a more significant presence of Radon gas in the sector [30,32].
The Machín volcano, with a crater 2.4 km diameter and 450m high, polygenetic, has a dome -developed through hundreds of thousands or million years-, with sedimentary and morphologic characteristics that suggest a Toba cone structure, built during the most recent phreatomagmatic [33]. According to studies related to the geological history of this volcano, a future pyroclastic eruptive episode would be deposited mainly in an area of 10 km 2 around the volcano edifice [34]. The similarity in the morphology of two volcanic structures and the activity of one of them must be taken into account, as one more element, in the evaluation of the danger of an eventual eruption. That could be the case with the similitude of the Cerro Machín Volcano and the Chichón or Chichonal Volcano in southern México in the state of Chiapas (17 o

Cerro Negro -Chiles
The Chiles and Cerro Negro are located on the Colombia-Ecuador border, 86 km from the city of San Juan de Pasto, at geographic coordinates 0 • 49'N 77 • 56'W and 0 • 46'N 77 • 57'W, respectively. This volcanic complex lies at the intersections of three faults: Chiles-Cerro Negro, Chiles-Cumbal and Cerro Negro-Nasta. The volcanic domes reach 4748 m a.s.l (Chiles), 4470 m a.s.l (Cerro Negro), and Figure 1: Artistic representation of 13 Colombian active volcanoes -Azufral, Cerro Negro, Chiles, Cumbal, Doña Juana, Galeras, Machín, Nevado del Huila, Nevado del Ruíz, Nevado Santa Isabel, Nevado del Tolima, Puracé, and Sotará-are displayed in three disperse clusters through the Cordillera Central. Because of their social significance and eruptive history, we shall briefly focus on four of them: Galeras, Nevado del Ruiz, Cerro Machín and Cerro Negro-Chiles. Galeras, Cerro Negro-Chiles volcanoes are found in the southern cluster, while Nevado del Ruiz and Cerro Machín are located within the most northern one. their craters have diameters of 1.0 km and 1.8 km, respectively. These two adjacent volcanic cones are collapsed towards the north (Chiles) and west (Cerro Negro), with the presence of geoforms of an already extinct glacier action. Their buildings are formed mainly by several episodes of lava and pyroclastic, with main volcanic products classified as andesites of two pyroxenes and olivines.
Although there are no historical records of eruptive activity, there is evidence of highly explosive stages, and the current activity is displayed by the presence of hot springs and solfataras. On the Ecuadorian side of Chiles, there is a seismological station which detects various activity.

Galeras
The Galeras volcanic complex -located in southwest Colombia: 1 • 13'18.58"N, 77 • 21'33.86"W-is the most active volcano in Colombia with the highest social risk due to its regular activity and the populated area that surrounds it.
Surpassing 5,000 years of antiquity, this volcanic complex has a base diameter of 20 km, a summit elevation of 4,276 m a.s.l., and a primary crater diameter of 320 m. The active cone, called Galeras Volcano, rises 1600 m above and approximately 9 km away from the city of San Juan de Pasto (capital of Nariño department) with a population of 313,000 inhabitants.
Galeras is characterized by andesite lava and pyroclastic with significant fallout deposits, displaying a conical shape with a large caldera at the top. After a long period of inactivity of more than 40 years, it awoke again in 1987, experimenting mainly minor eruptions, some with explosive character: fumarolic formation and enlargement, strong tremors, shockwaves and emission of pyroclasts and ashes [35,36]. Since 2009, the activity has been considerably reduced to the expulsion of ashes -columns that have reached 10 km-and shockwaves [37].

Nevado del Ruíz
The ice-covered volcano Nevado del Ruíz (NRV), found at the Cordillera Central -4 • 53'N and 75 • 19'Whas an altitude of approximately 5390 m. a.s.l., covering an area of more than 200 square kilometers. Its main crater (Arenas) is one kilometer diameter and 240 m deep. La Piraña and La Olleta are two small parasitic edifices, and there are four U-shaped amphitheaters produced by flank collapse and fault activity [38].
This enormous volcano is located at the junction of two fault systems: the N75 • W normal Villa María-Termales fault system and the N20 • E right-lateral strike-slip Palestina fault system [39]. The north and northwest borders show uneven geometries caused by the location of large amphitheaters on the upper part of the volcano, the southern and southwestern fringes are marked by sharp regular slopes while the east and southeastern fringes present moderate-to-strong declivities and a significant thickness of glacier deposits.
The eruptive history of the Nevado del Ruíz runs from the Pleistocene to the present, and its stratigraphy has three main stages related to the alternate construction-destruction of its edifice: Ancestral Ruiz (2-1 mega-years), Older Ruiz (0.8-0.2 mega-years) and Present Ruiz (<0.15 mega-years) [40]. The present emplacement of lava domes is made of andesite and dacite inside older calderas [41]. During the past 11,000 years, this volcano has passed through at least 12 eruption stages, which included multiple slope failures (rock avalanches), pyroclastic and lahars, leading to the partial destruction of the summit domes [40,41]. The eruptions of the last thousand years have mostly been small, excluding some like the phreatic-magmatic eruption on November 13, 1985, which involved the partial melting of the glacier cap and consequent lahars, which reached and destroyed the municipality of Armero-Tolima and caused a large number of casualties.

Muon flux simulation chain
Particles measured at the ground level (secondaries from now on) come from a chain of interactions and reactions, started by the primaries impinging upon the outer atmosphere. The modulation of secondaries needs to be monitored and carefully corrected by taking into account atmospheric factors that could modify its flux. Thus, a complete and detailed simulation chain -considering factors such as geomagnetic conditions, atmospheric reaction and detector response-is needed to characterize the expected flux at the detector level. Any attempt to estimate the expected flux of secondaries at the detector level should consider a detailed simulation that takes into account all possible sources of flux variations of processes occurring at different spatial and time scales. We can illustrate this conceptual scheme as [42,43,44]: →Atmosphere Secondaries → Detector response → Signals.
We start our simulation chain by characterizing the effects of the geomagnetic field (GMF) on the propagation of charged particles contributing to the background radiation at ground level. This is included through the calculation of directional rigidity cut-off, R c , at the detector site, determined by using the Magnetocosmics code, which implements the backtracking technique [45,46]. GMF at any point of Earth is calculated by using the International Geomagnetic Field Reference [47], for modeling the near-Earth magnetic field (r 5R ⊕ ) and by the Tsyganenko Magnetic Field model version 2001 to describe the outer GMF (r 5R ⊕ ) [48].
The second link in this chain corresponds to a detailed simulation of extensive air showers produced during the interaction of the flux of primaries -from protons to irons, i.e., 1 ≤ Z ≤ 26-through the atmosphere. To obtain this very comprehensive secondary flux at ground level we use the Corsika code [49], with a specially tuned set of simulation parameters adapted to a particular geographical site. We will display these parameters in Section 5.1 for the particular case of the Machín volcano. For the hadronic model of interaction in the distribution of secondaries at the level of the detector, we use QGSJET-II-04 [50] as a model, and for low energies, we choose the default option, GHEISHA-2002 [51].
To identify possible muon observation sites in Colombia, we shall implement the previously described calculation scheme only at secular geomagnetic conditions -i.e., static geomagnetic corrections-focusing on the detailed calculation of the crossing muon flux and the stopping power of the volcano edifice at different volcano sites. Although we have implemented this general scheme, we found that geomagnetic corrections seem not to be significant for the geographical location of Colombian volcanoes, because it mainly affects low energy primaries [44].
Finally, the detector response to the different types of secondary particles at ground level is simulated using a Geant4 model [52,53] for both the scintillator panels and the water Cherenkov detector, but this last step of the simulation chain will not be considered here but will be detailed in a coming work [54,55,56,57].

Directional muon rock opacity
The open sky flux estimations at the ground level, influenced by various environmental parameters -altitude, geomagnetic corrections, solar modulation, and atmospheric variations-induce some critical features in the instrument design.
A comparison of the open sky particle flux Φ os with the flux, Φ, emerging after crossing the target provides information about its density gradient. To estimate this target density gradient we define the directional rock opacity as the mass density ρ integrated along the muon path L as: where ξ is a characteristic longitudinal coordinate through the volcano, L is the total distance traveled by muons in the rock,ρ is the average density within the volcano. In our approach, we are assuming that the trajectories of muons are straight lines which are not affected by Coulomb scattering processes; therefore density distributions within volcano edifices can be inferred from the variation of its flux, ∆ = Φ − Φ os . A more detailed calculation, including second-order effects, is being carried out and will be included in future characterizations of the selected places.
The muon energy loss along each path can be modeled as for each muon arrival direction considering a uniform density distribution. Here E is the muon energy; a(E) and b(E) are functional parameters depending on the rock composition/properties and (L) is the density integrated along the trajectory of the muons (the opacity defined by eq. (1)). The coefficient a(E) represents the energy loss due to ionization, while b(E) takes into account the contribution of radiative losses, mainly Bremßtrahlung, nuclear interactions, and pair production. The main parameters to estimate the coefficients a(E) and b(E) are the average ratio < Z/A > between the atomic and mass numbers of the material [58,59].

MuTe: Colombian Muon Telescope
There are three main types of detectors implemented for volcano muography: (a) nuclear emulsion detectors, (b) scintillation detectors, and (c) gaseous detectors. Each one has its pros and cons as described in [8,12]. The main issue in muography of large-size objects is the flux overestimation due to the background sources such as the soft component of Extended Air Showers (EAS) [60,61] and upward-going muons [10]. Here, we shall present a new hybrid Muon Telescope (MuTe), combining the facilities of a hodoscope with a water Cherenkov detector (WCD) in order to solve that drawback. The hodoscope estimates the event flux per trajectory depending on the fired pixels in the front and the rear panel. A Time-of-Flight, ToF, system measures the time taken by the crossing particles subtracting the nano-scale time stamp recorded in both panels. The WCD senses the Cherenkov photons generated in the water due to the interacting charged particle. The recorded photon yield is equivalent for the energy loss.

The instrument
Just for completeness, we briefly describe here the muon telescope its acceptance and noise reduction capabilities, later we employ these characteristics to estimate the time exposure of the instrument at the observational sites. A detailed description of the instrument capabilities will be discussed shortly elsewhere [56,57]. Figure 2 illustrates MuTe design which combines two detection techniques: a hodoscope formed by two detection planes of plastic scintillator bars, and a WCD, in an innovative setup which differentiates it from some previous detectors.
• Water Cherenkov Detector: The WCD is a purified water cube of 120 cm side, located behind the rear scintillator panel (again, see figure 2), which acts as an absorbing element and as a third active coincidence detector is capable of isolating the muonic component of the incident particle flux. From the charge histogram, obtained by time integration of the individual pulses measured in the WCD, it is possible to separate two components of the incident flux: electromagnetic part (photon, electron & positron) and the µ−component [42]. Additionally, due to its dimensions and location, it filters most of the electrons backward noise while protons/muons moving backward -which could cause overestimation in the hodoscope counts [64]-are rejected by a Time-of-flight system.
Muon generated events deposited in an energy range of (144MeV< dE/dx < 400MeV), represent only about 40% of the WCD-hodoscope acquired events. The other 60% of data is composed by (e ± ) events under 144MeV and multiparticle events above 400 MeV. Subsequently, low-momentum muons (< 1 GeV/c), which are scattered by the volcano surface, are also rejected.
The Colombian MuTe combines particle identification techniques to discriminate noise background from data. It filters the primary noise sources for muography, i.e., the soft-component (e ± ) of Extensive Air Showers (EAS) and scattered/upward-coming muons. Particle deposited energy identifies Electrons/positrons events in the WCD, while scattered and backward muons are rejected using a pico-second Time-of-Flight system.

Telescope acceptance
The acceptance of the instrument is a convolution of the telescope geometry (number of pixels, pixel size, and panel separation) and it is obtained multiplying the detection area by the angular resolution, i.e. T (r mn ) = S(r mn ) × δΩ(r mn ), where r mn represents each discrete muon incoming direction, which for an array of two panels with N x × N y pixels we can identify (2N x − 1)(2N y − 1) different particle trajectories [17]. Moreover, the number of incident particles N ( ) can be expressed as and I( ) is the integrated flux (measured in cm −2 sr −1 s −1 ), T represents the acceptance (measured in cm 2 sr) and ∆T designates the time exposure (in seconds). It is possible to obtain a simple relationship between the exposure time and the desired opacity resolution as with ∆I( 0 , δ ) the flux variation due to the different opacities 0 and 0 + δ , and c is a parameter measuring the confidence level in terms of the standard deviation of the measurement. The above expression gives a bound for the minimum exposure time needed to distinguish opacity differences across the geological object [17].
In figure 3 the angular resolution and the acceptance function for our telescope (with 900 pixels and 3481 discrete r mn directions) is shown. The total angular aperture of our telescope is roughly 582 mrad with the maximum point at 1.6 × 10 −3 sr and, as expected, the largest detection surface corresponds to the normal direction r 00 , reaching ≈ 6 cm 2 sr.

Machín, Chiles, Cerro Negro and Galeras Volcanoes
In this Section we shall first present a detailed study of the muon propagation across the Machín volcano and calculate the corresponding exposure times for four observation points. We shall also display estimations of the outgoing muon flux from the Machín dome that could be detected by our muon telescope. Secondly, we also display a ray-tracing analysis for muon propagation for three other active Colombia volcanoes: Chiles, Cerro Negro, and Galeras.
From the information gathered, we identify some critical parameters that could limit the application of muography in Colombian inland volcanoes and devise a "rule of thumb" criteria to identify possible muography volcano candidates.

Detailed study of the flux from the Machín volcano
We start our study with Machín volcano because it is one of the most dangerous active volcanoes in Colombia. It has an average height of 2750m.a.s.l; a crater of 2.4km of diameter blended into the landscape of its nearby topography, which makes it practically invisible. This fact increases the risk of the Machín volcano for the surrounding population. Additionally, the morphological similarities between Chichonal and Machín volcanoes -both are considered stratovolcanoes and went unnoticed as volcanoes for decades [65,66]-is impressive. The resemblances of these two volcanic buildings amaze and, the chipping of the eruption of the Chichonal volcano -the largest in the history of México, with an affectation of almost 100km around-should alert about the potential danger of Machín.

Cosmic ray spectral composition at the Machín volcano
We implement the simulation chain described in section 3.

Muon propagation through the Machín volcano
Next, we calculate the differential flux of muons at the above mentioned four observational points as a function of the direction of arrival and determine the maximum possible depth that can be observed from each point. This is shown in the left plate of figure 5, where we display the momentum spectra for muons emerging in five angular bins after crossing a standard variable rock with density 2.65 g/cm 3 , measured at the point P 1M of Machín volcano. As it can be clearly appreciated, the muon flux decreases considerately reaching 1 cm −2 sr −1 day −1 for p µ 500 GeV/c.
We have employed the topography from the NASA Shuttle Radar Topography Mission global digital elevation model of the Earth 2 , with resolution of 90m×90m (see left plate in figure 6 for the particular case of the Machín volcano). Next, we calculate all the possible distances crossed by each muon path integrating equation (2) and (1) for standard-rock-model dome, with the coefficients a(E) and b(E) obtained from the Particle Data Group [58] 3 .
At Cerro Machín we have identified four observation points: P 1M , P 2M , P 3M and P 4M around it (see Table 1). Figure 6 displays the ray-tracing technique implemented for those points with the corresponding muon propagation distance through the topography, as well as the angular distribution of these distances around the upper part of the volcano. Emerging muons with these trajectories have crossed about 600 meters of rock within the geological structure.  We have found that very few muons -the most penetrant component of the particle shower, ranging from tenths to few thousands of GeV/c-, ≈ 10 −1 muons per square meter, per day, (see figure 5, left plate)-can cross almost 1, 000m of standard rock (see figure 5, right plate or figures 6 ). Therefore we set 1, 500m as the upper bound for distances that can be traveled by the most energetic horizontal (Zenith angles 70 degrees) muons at any Colombian Volcano. Our analysis concludes that starting from a few GeV/c there is no significant effect of the geomagnetic correction on the muon flux at any geographical zone in Colombia, but this correction is, in general, essential to determine all the possible particle background flux at other sites with different latitudes [44].
As we have mentioned above, we set a minimum flux of 100 µ/pixel and by using equation (4) we estimate the minimum time exposure needed to examine the inner structures of the volcano edifice.
In figure 8 we sketch contour lines representing the exposure times required, depending on the point of observation. If we assume an acceptance of 6 cm −2 sr., we will need at least 100 days (∼ three months) to explore a depth of 150 − 160 meters.
As explained before, exposure times, opacity (directional average density) and instrument resolution are related through equation (5) [17]. This is another variable: the expected exposure times needed to resolve average density differences of ≈ 10%, are shown in figure 9, for the zenithal range 66 • < θ < 84 • .

Ray-tracing for Chiles, Cerro Negro and Galeras
The idea in this Section is to look for potential observation points where the muon flux travel distance, inside the geological structures, is less than 1, 500m. It should be emphasized that almost all Colombian inland volcanoes are surrounded by other geological structures that screen the atmospheric muon flux.
Thus only a few potential observation points are available.
We examine the other three critical Colombian volcanoes using the ray-tracing technique based on the topography, already mentioned, above of the NASA global digital elevation model. First, the code returns the latitude and longitude coordinates in vectors within a matrix of elevation values and redefines the reference system (decimal degrees) to a local reference system (meters). Then we trace the muon trajectory considering in detail the topography around each volcano.

Observation points at Chiles and Cerro negro volcanoes
At Chiles we have identified four points (see Table 2) and in figure 10 we can see the ray-tracing technique implemented for P 1Ch , P 2Ch , P 3Ch and P 4Ch to determine the distances of muon propagation through the volcano, as well as the angular distribution of these rays around the upper part of the volcano. Although the distances traveled by the muons to the four observation points are less than 1, 500m, the volcano was discarded due to the difficult access to these points.  In the case of Cerro Negro we have identified four points around it (see Table 3). In figure 11 we  can see the ray-tracing technique implemented for P 1CN , P 2CN , P 3CN and P 4CN to determine the distances of muon propagation through the Cerro Negro volcano, as well as the angular distribution  : Exposure times for observation point P 1M at Cerro Machín needed to identify differences of −10% in the averaged directional density for different zenith angles and telescope acceptance. We obtained exposure time lapses between two days and up to more than six months to achieve the desired density resolution, at different zenith angles.
of these rays around the upper part of the volcano. Notice that to obtain a reasonable muon flux the telescope should be tilted about seven degrees making the depth of investigation very short.
Although the distances traveled by the muons measured from points P 1CN -P 4CN are less than 1, 500m, this volcano, as well as Chiles, are discarded due to the difficult access to their potential observation points. The Chiles and Cerro Negro volcanoes are on the Colombia-Ecuador border, an intricate zone, where the safety of the scientific personnel and equipment is not guaranteed by either country. Figure 10: Particle trajectories crossing the Chiles volcano structure to observation point P 1Ch , P 2Ch , P 3Ch and P 4Ch . Points P 1Ch and P 2Ch are in areas that are difficult to access and points P 3Ch and P 4Ch in high risk areas Figure 11: Particle trajectories crossing the Cerro Negro volcanic structure to the observation point P 1CN , P 2CN , P 3CN and P 4CN .

Observation points at Galeras volcano
Galeras volcano has been active for more than a million years. Its most recent phase of activity began about 4500 years ago and included six significant eruptions. After quiescence dating from 1948, Galeras renewed its activity in 1988 with the emplacement of a dome, which was destroyed in an explosion in 1992 [69]. On January 14, 1993, an explosion in Galeras crater killed six visiting scientists and three tourist [70].
The Galera's intense activity and its proximity to an important city, generate a significant interest to study this Volcano with all techniques available and muography do not escape from this pressure. Several local researchers have made a notable effort in modeling muon transport properties across this dangerous volcano [23,27,28], and some of them even proposed a design for a muon telescope [27]. It is not clear which are the geographical coordinates where they intend to place the instrument and what will be the typical exposures times so as to have statistically reasonable counts.
We have proceeded to identify where safe observation points might be to install our MuTe, and how long would it takes to obtain a valuable density distribution. Following this procedure, we have identified four points around Galera Volcano (see Table 4 and figure 12). However, all are in the high volcanic risk area, and so located that the detected muons through the amphitheater and crater overlap, making difficult the analysis (see figure 13.

Volcano observation site determination criteria
Based on the detailed analysis of the Machín volcano, and the ray-tracing studies of the Cerro Negro-Chiles complex and Galeras, we have devised a "rule of thumb" criteria to select the tentative muography observation sites of Colombian volcanoes.
These criteria for determination of the muon observation sites for the Colombian mainland volcanoes, -surrounded by other geological structures that could screen the atmospheric muon flux and having insecure access to various regions because of the internal conflict-is qualitatively different from those made on peaceful islands with volcanoes free from the screening of other mountain systems. Thus, to determine muon observation points for active volcanoes in Colombia, we have devised a mix of technical and logistic criteria, which we call the "rule of thumb" criteria that should be fulfilled by the potential sites and which are listed below: • Criterion 1: At the observational level, is the volcano base width less than 1, 500 m? This criterion is necessary because the energy of atmospheric horizontal muons are two orders of magnitude lower than vertical muons and, these incident particles need two more orders of magnitude to cross the volcanic structure [64]. Therefore, most energetic horizontal muons can only cross 1, 500m of standard rock, and the horizontal line of sight from the tentative observation point should be less than this distance. This criterion can be appreciated quantitatively for trajectories/flux in the case of Machín volcano (figures 6 and 7), where only few (the most energetic and horizontal) muons can cross the structure while traveling distances greater than 1, 500m.
• Criterion 2: Are there tentative observation points where the surrounding topography does not interfere with the target? Muons impacting the telescope should cross only the structure under study. Nearby mountains and any other geological formations neighboring the target volcano, must not contribute to the opacity. This requisite imposes a severe restriction on the tentative view points for the few observational places where a small window is present, with no mountains or other screening geological structures.
• Criterion 3: Are the sites accessible and secure? Sites must be easily accessible, and the telescope should be securely transported and placed on the field. It is essential to consider: the weight and size of the assembled telescope and its parts; also the quality and accessibility of water resources in the area. Additionally, the volcano to be studied should not be cataloged in a situation of abundant activity due to the danger of volcanic products and processes associated with eruptions such as ashfall, pyroclastic materials, lahars, floods, among other risks, as well as earthquakes and landslides that may cause severe damage to instrument and personnel. Last but not least the -vanishing-internal conflict persisting in several regions of the Colombia countryside and safe access should be taken into account.
We devised these criteria based on the detailed study of Machín volcano and the ray-tracing analysis for 13 Colombian volcanoes. The results of the application of the above criteria are summarized in table 5 and lead to the conclusion that the only Colombian volcano that could be studied through muography is Cerro Machín [71]. We hope that shortly with the end of the internal conflict, Cerro Negro and Chiles will be safely accessible because they are good candidates to be studied with muography, but today they are not yet freely reachable.

Discussion and Conclusions
In this work, we present the first comprehensive simulation of muon flux for the Machín volcano dome.
We have also carried out a ray-tracing analysis for 12 other inland volcanoes in Colombia surrounded by complex topographic environments. After a detailed study from the topography, we have identified the best volcano to be studied, spotted the finest points to place a muon telescope and estimated its time exposures for a significant statistics of muon flux.
The rationale of our new approach stems from a four-step methodology: 1. A "rule of thumb" criteria. We have devised a mix of technical and logistic criteria -the "rule of thumb" criteria-and applied them to 13 Colombian volcanoes. We have found that only Cerro Machín, located at the Cordillera Central (4 • 29'N 75 • 22'W), can be feasibly studied today through muography. Cerro Negro and Chiles could be good candidates shortly.

2.
A unabridged simulation of open sky particle spectrum and composition. Corsika has calculated the energy spectrum and composition of open sky secondaries at Cerro Machín and filtered them within Magnetocosmics framework, providing a detailed description of different types of particles in the MeVs to TeVs secondary energy range [42,43].

3.
A detailed calculation of the emerging muon flux.
With the above open sky particle flux and precise topographical information surrounding the volcano, we have simulated the muon propagation inside the geological edifice and estimated the emerging muon flux at four different points around Cerro Machín. This was carried out integrating the energy loss equation 2, with the coefficients a(E) (ionization) and b(E) (radiative losses) from reference [58].
4. An estimation of the telescope time exposure. With the emerging muon flux, we have calculated the time exposures of the instrument for an arbitrary statistics of 100 events at each pixel and, different values of the telescope acceptance (see Figure 9).
As we have mentioned before, most of the previous muography studied volcanoes -Mount Asama [18] in Japan; Puy de Dôme [72] in France; Mount Etna [9] in Italy; La Soufuriere [73,74] in Guadalupe, to mention the most relevant studies-are topographically isolated with a relatively good and accessible observation points. None are surrounded by geological structures screening the scarce high energy horizontal muons. However, Colombia -and surely all other Andean volcanoes-is very different; most of the active volcanoes are along the Cordillera Central, surrounded by higher altitudes shielding cosmic ray flux. Therefore, we developed a methodology to identify the most feasible candidates, and only Cerro Machín emerged as a possibility.
Instead of using phenomenological and pseudo-empirical formulas to estimate the background flux at the volcano site (see [18] and references therein), we proceeded to simulate its spectrum and composition, at each particular geographical location, with two standard astroparticle tools: Corsika and Magnetocosmics. We found that incident muons range from 0.1 GeV/c to 10 TeV, and the flux of high energy muons is very feeble:≈ 10 muons per square meter per day at zenith angles θ ≈ 82 • − 84 • .
With the above simulation as an input, and including precise topographical information, we calculate the propagation of muons through the geological edifice and determine the emerging muon flux that could be detected at several particular observation points around the volcano (see figure 6).
We have simulated more than 50 days of muon flux to estimate the minimum time needed to obtain a flux of 100 muons per pixel and found out that we require at least three months of data acquisition to explore the Machín dome (depending on the observation point) at a depth of ≈ 110m.
Then, to discriminate density variations of 10%, we evaluated time exposures for our hybrid instrument as function of the acceptance. With these preliminary simulation results and by considering the standard configuration of our telescope, we have estimated time intervals between 100 to 125 days for the upper end (114 − 150 m) of the volcano edifice. These results were recently reconfirmed using more precise muon underground propagation codes [67,68].
Muography can not image deep volcano structures, but it seems to be useful in determining shallow phenomena with an excellent spatial resolution. This technique can not give direct information on when a volcano will erupt, but it could provide significant insights about possible eruption processes, in the upper end of the edifice. This emerging technique requires significant progress in data analysis, treatment and interpretation of the experimental data obtained. For a bright future, it depends on the synergy between two active international communities: particle physicists and geophysicists [12].

Acknowledgments
The authors acknowledge the financial support of Departamento Our simulation codes can be found at https://github.com/AstroparticulasBucaramanga and other pertinent data can be also obtained from https://zenodo.org/record/807741#.WUMdlMaZPex