Research Article  Open Access
Planning Tunnel Construction Using Markov Chain Monte Carlo (MCMC)
Abstract
Tunnels, drifts, drives, and other types of underground excavation are very common in mining as well as in the construction of roads, railways, dams, and other civil engineering projects. Planning is essential to the success of tunnel excavation, and construction time is one of the most important factors to be taken into account. This paper proposes a simulation algorithm based on a stochastic numerical method, the Markov chain Monte Carlo method, that can provide the best estimate of the opening excavation times for the classic method of drilling and blasting. Taking account of technical considerations that affect the tunnel excavation cycle, the simulation is developed through a computational algorithm. Using the Markov chain Monte Carlo method, the unit operations involved in the underground excavation cycle are identified and assigned probability distributions that, with random number input, make it possible to simulate the total excavation time. The results obtained with this method are compared with a real case of tunneling excavation. By incorporating variability in the planning, it is possible to determine with greater certainty the ranges over which the execution times of the unit operations fluctuate. In addition, the financial risks associated with planning errors can be reduced and the exploitation of resources maximized.
1. Introduction
Underground mining represents a fundamental pillar of ore production in Chile. It is assumed that in the coming years the proportion of underground mining compared with opencast mining will increase as mineral resources accessible to surface exploitation become progressively exhausted.
One of the main activities involved in underground mining is tunnel construction, or, more generally, horizontal works, because this produces the infrastructure that provides access to the ore for extraction. Here, a “tunnel” should be understood as any underground excavation whose purpose is to join two points.
Given the importance of tunnels for mining, it is evident that there is a need to have a methodology that allows accurate planning of their excavation. To achieve this goal, the Markov chain Monte Carlo (MCMC) method is appropriate, since the construction of a tunnel is a cycle of activities consisting of unit operations, each of which exhibits a variability that can be represented in terms of a probability distribution function (PDF). Furthermore, the success or failure of the construction cycle is related to its actual duration compared with what was planned, which also depends on the time at which the cycle begins within the day’s work shift. Thus, the construction cycle of a tunnel is dependent on the success or failure of the immediately preceding cycle, and therefore the event’s probability of success is related to its predecessor, constituting an MCMC relation [1].
2. Tunnel Construction
There are several methods of tunnel excavation. This paper will focus on the construction of tunnels by drilling and blasting [2]. This technique involves an excavation or work cycle comprising a number of different activities. Suorineni et al. [3] mention the following unit operations: drilling of the tunnel surface, loading of explosives and blasting, ventilation (considered as an interference within the cycle), scaling and loading of the blasted material, and fortification (bolts, nets, and shotcrete, among others). Figure 1 illustrates the drilling and blasting cycles involved in tunneling.
The aim of the excavation cycle is to break up the rock with explosives, giving the required crosssectional shape while the tunnel advances proportionally to the length of the drilling in the tunnel face. In this way, with successive cycles, the infrastructure is built gradually until the tunnel has been completed. However, even with knowledge of the number of unit operations in each cycle, and the time generally taken to perform each one, it is very difficult to determine exactly the total time required to complete the construction of the tunnel, mainly because all of the operations are subject to variations that depend on unforeseen events (although they can be associated with a probability of occurrence).
On the other hand, the work cycles, and therefore each of the unit operations, are executed within welldefined time periods (work shifts), which are framed within a 24hour period. Usually, mining operates in continuous time periods without stoppages in production, and therefore tunnel construction proceeds in the same continuous manner. This is particularly important because the relation between the duration of the construction cycle and the period defined for the work shift will affect the efficiency of the cycle. Some of these inefficiencies result from changes in work shift that interfere with the working cycle. That is why Chilean law [4] specifies several types of work shifts, with various configurations as shown in Table 1. The choice among these is made on the basis of the estimated duration of the tunnel’s construction cycle. For example, if the cycle time is estimated as less than 8 hours, the work shift that fits this best should be used, in this case T1 (Table 1), because this allows three cycles per day and thus a more rapid advance of the tunnel.

3. Planning of Tunnel Construction
Currently, to plan the construction of a tunnel, whatever its purpose, fixed values of the relevant parameters are used, giving consistent results. This is reflected in the following equation, which gives the construction speed of the tunnel in days related to the drilling length unit () in terms of the drilling length (), weighted by the effectiveness of the blasting (), divided by the sum of the times for the unit operations in hours (), and divided in turn by a factor that involves the unproductive times () in relation to the 24 hours of the day: Then, with knowledge of the length of the tunnel (), the execution time () is given byThe results obtained with this approach have until now always been used to plan this type of construction, but they can be improved by including the variability of each of the unit operations involved.
4. Monte Carlo Method and Tunnel Construction Planning
As already mentioned, tunnel construction involves excavation cycles consisting of unit operations that can be represented by PDFs, and it is clear that this process can be simulated using a Monte Carlo method [5]. With this approach, an excavation cycle () is simulated as follows:where is the inverse probability function of each of the unit operations and is a pseudorandom number between 0 and 1.
Assuming that the number of excavation cycles required to construct the total length of the tunnel is known, since the advance achieved in each cycle is determined by the drilling length, which remains unchanged throughout the construction, the theoretical time taken for construction () is given by the following equation:where is the number of cycles. For the model presented in (4) to represent reality, it is necessary to include in the construction of the PDFs the unproductive times and inefficiencies associated with each activity.
However, the duration of tunnel construction cannot be estimated using the Monte Carlo method alone, because this method does not take account of an aspect that is extremely important, namely, the fact that the probability of success of a cycle depends on the preceding cycle.
5. Application of Markov Chains
As already mentioned, the use of Markov chain theory is appropriate in this context, considering the characteristics of tunnel construction. Construction in mining takes place continuously 24 hours per day, and in general the working day is broken up into two or three periods (shifts), depending on the chosen workday (see Table 1) and as permitted by Chilean law [4]. Taken into account this work structure, tunnel construction in mining is faced with some particular problems that are mainly the result of inefficiencies due to changes of work teams and their transfer to the working areas.
In underground mining, owing to the specific characteristics of the work, excavation cycles are generally kept as multiples of working shifts, for a duration of less than 8 hours, for example, with preference being given to a T1 type of shift over T2 or T3 (Table 1), so that three cycles per day can be run.
It is possible that, owing to particular aspects of operational interference or inefficiency, an excavation cycle will not fit the established workday, increasing the duration of the cycle and affecting the next cycle. On the other hand, if success in the execution of an excavation cycle is represented by its completion within the established work shift or group of work shifts (with the cycle otherwise being considered a failure), then this condition in turn reduces the possibility of success of the following shift, because it takes time away from the latter and furthermore adds unproductive time due to activities interrupted as a result of the change of work. These situations can be considered as processes that can be modeled using Markov chain theory [1].
The algorithm presented in (4) cannot model such behavior, because it is unable to determine the simulation time of an excavation cycle as a function of the duration of the work shift. Therefore, to incorporate this behavior, it is necessary to have an algorithm for evaluating this simulation time.
6. Simulation Algorithm
To predict tunnel construction time, the Monte Carlo method appears to be an appropriate tool to use together with the Markov chain principle, given that it is a stochastic simulation that allows analysis of complex systems with several degrees of freedom. The Monte Carlo method [5] has become one of the most common ways to solve complex mathematical problems by random sampling [6–10]. It consists in generating random or pseudorandom numbers that are entered into an inverse distribution function, delivering as a result as many scenarios as the number of simulations performed [11]. The estimation will be the more precise the greater the number of iterations that can be done.
To use the Monte Carlo method, the unit operations are identified and each is assigned a PDF that depends on its nature and on the results of field sampling.
If the inverse functions of the PDFs of each unit operation of the excavation cycle are fed with random numbers, they will give as a result the duration of each operation. If the times thus obtained are added, this gives the total duration of the excavation cycle.
Once we know the duration of the excavation cycle, we must also consider another very important variable, namely, the distance advanced, or the real advance, after blasting (). This distance can also be described by a PDF, since it corresponds to the drilling distance () as affected by the blasting efficiency (), as shown in (1). The drilling distance is a fixed value that depends on the characteristics of the drilling equipment, but the efficiency of the blast depends on the condition of the rock, variations in geological structure, and the characteristics of the explosive used, among other things, making this parameter vary from one blast to another.
Thus, if the durations of all the excavation cycles and the corresponding distances advanced are known, it is possible to determine the time taken for tunnel construction. Simulating this as many times as possible, a large number of scenarios are produced, generating a PDF of the time taken for tunnel construction.
The algorithm (see Algorithm 1) is composed of three loops, which control the number of simulations required, the required tunnel length, and the existing relation between the duration of the work shift and that of the tunnel excavation cycle. This point is fundamental for this work, being very important when it comes to choosing the best shift configuration to use. All these items are necessary to allow the simulation of the total construction time.

The proposed scheme consists of three inclusive loops dependent on each other. The operating form is that the first loop, which contains the other two, controls the number of required simulations, on the basis that each simulation is the construction of a tunnel with specified length.
The second loop imposes the condition that the construction does not exceed the defined tunnel length, with every advance being estimated by the PDF of the performance of the blast multiplied by the drilling length. The drilling length is a fixed value, and is added consecutively until the required tunnel length is achieved.
Finally, the third loop has the function of adding consecutively the times for the unit operations in each cycle and comparing this total time with the established work shift, a fundamental aspect of this work.
This last loop is the key to the simulation, because it constructs the cycle within the shift. This procedure is carried out using the Monte Carlo method, where the inverses of the PDFs associated with the execution times of the unit operations are applied. The times obtained from this simulation are added, and it is determined whether the cycle can finish during the operating shift.
The proposed simulation algorithm is detailed in the following section.
6.1. Control of the Number of Simulations
As already mentioned, the function of the first loop is to control the number of simulations required, taking into account that each simulation estimates the time required to build a tunnel with an already defined length. This loop, which contains the other two, starts with the variable “,” which represents the number of the simulation under way and is initially assigned the value 1. The first loop is “WHILE ,” where “nsim” corresponds to the number of required simulations.
Some variables are then defined to start the execution of the algorithm. The variable “dev” is an auxiliary adder that is assigned an initial value of 0 and is increased in relation to the advances per blasting at the end of the second loop. The purpose of this variable is to control the simulation in relation to the tunnel length built, executing the second loop until the desired length “WHILE ” is achieved, where “ltunnel” represents the length of the tunnel to be built.
In the first loop, there is the variable “operation,” whose purpose is to identify the operations that will be performed in the third loop and that are adapted according to the operations of the construction cycle described by Suorineni et al. [3]. It should be noted that this is the only alphanumeric variable. The value assigned by default in this location is “op.1,” since it indicates the first operation of the cycle, which is represented by the suffix 1. Every time that we begin simulating the construction of the tunnel, we will start with drilling, the first operation of the cycle.
On the other hand, the variable “shift” defines the work shifts required for the construction of the tunnel, and it is a counter that is modified at the end of the second loop, storing the data in a matrix “” that is of dimension equal to the number of required simulations and that will provide the data for later analysis.
Finally, the variable “,” an auxiliary variable used for storing the sum of the times of the operations required for the tunnel construction, is set equal to zero every time the tunnel construction starts, but it can also be initialized, depending on conditions that will be explained later, within the second loop.
6.2. Control of the Simulated Construction Advance
The second loop imposes the condition that the simulated construction advance does not exceed the proposed length, and to that end the construction takes place in the third loop, whose purpose is to control the duration of the cycle in terms of the duration of the work shift “WHILE .”
The verification expression is true as long as the sum of the times of the cycle’s operations “” is less than the duration of the shift “” minus a tolerance time “tol.” This last variable is an operational parameter that indicates if it is possible to continue with the next unit operation within the shift or if the operation is to be passed to the following shift.
6.3. Calculation of the Time for Each of the Unit Operations
The cycle’s duration in the work shift is built successively in a selection routine “CASE” that adds the time for each operation until the end of the cycle; then, the following cycle can start again within the same shift or stop the execution to retake it on the following shift, and this depends on the operational tolerance “tol” that is estimated for the execution of the tunnel. We will go into this point more deeply when we apply the algorithm.
The “CASE” routine in the third loop has the function of arranging the operations so they take place one after the other and also of evaluating the time taken under the loop’s condition, in order to see if this is still within the duration of the chosen shift.
The times for each operation belong to the probability distributions used to represent the process. In our case, the variable “DI.op.n” corresponds to the inverse distribution of the operation specified in the suffix, in this case “,” and this variable is a function of random numbers between 0 and 1.
By means of the variable “rand#,” which represents random numbers between 0 and 1, the values of the operation are generated and fitted to the distribution used, as pointed out by Sobol [11]. Once the operation has been executed, the variable “operation” stores the value of the following operation so that, in the next iteration, as a result of “CASE,” it keeps advancing.
At the end of the third loop, there is a conditioning routine depending on whether the cycle ends together with the shift or is interrupted. This routine evaluates whether “” is greater than “,” telling the algorithm whether the next shift should add an activity restarting time “,” where “beg” corresponds to the restarting time, which is not included in any of the unit operations. If the shift ends cutting an activity, this restarting time is added. In the opposite case, where the activity ends within the shift, it is not necessary to add the restarting time. If appropriate, this restarting time will be added to the next sequence of the loop in the corresponding operation.
Also, at the end of the second loop, a work shift is added in the variable “shift,” and the advance is added in the variable “dev” only if it has gone through the last operation where the advance caused by the blasting, “,” is found, where “av” reflects the advance of the tunnel (in meters) and “DI.rec(r#)” is the inverse distribution of the percentage efficiency of the advance caused by the blast.
In this way, the successive simulations are constructed, delivering the time taken for each simulated tunnel, and accounted for in work shifts.
Taking in account the possibility of iterating as many times as necessary, we will have a representative sample of the population from which we can infer the most probable duration of the tunnel’s construction.
7. Application of the Algorithm to Tunnel Construction
Minera San Pedro Limitada (MSP) has several copper ore deposits in the Lohpan Alto district, located in the Coastal Range of Central Chile. One of these deposits is Mina Romero, where the ore will be removed by underground mining [12].
To gain access and prepare the mineralized body for its exploitation, MSP has planned the construction of a 560 m access tunnel with no slope and in a straight line, with a cross section of approximately 3.5 m × 3.0 m.
The equipment used for drilling is an electrohydraulic drill of 45 mm diameter and the explosives are ammonium nitratefuel oil (ANFO) and Tronex (a derivative of dynamite), initiated by a nonelectrical shock tube detonator (NONEL). The mucking equipment is a load haul dump (LHD) of 6 yd^{3}.
To estimate the duration of construction of this tunnel, MSP, based on its experience in similar projects, has considered that, for every three work shifts with a duration of 9 actual hours each, it is able to carry out four cycles with a drilling length of 1.8 m, giving a rate of advance of 2.4 m per shift, so, with two work shifts per day, an advance of 4.8 m/day would be achieved.
Taking the above figures into account, MSP has estimated that the project should take 117 days, although previous experience in mines close to Mina Romero has shown that this estimate is not precise, because there are often delays that have not been considered at the time of planning.
To apply the proposed methodology, data from areas with similar geological and operational characteristics to those that will be faced in Mina Romero have been used. For that purpose, exploratory tunnels have been made in the upper part of the deposit, with the same cross section of the one that will be built and in rocks with similar characteristics to those of the Mina Romero access tunnel.
The cycle has been divided into five activities: drilling, loading and blasting, ventilation, scaling, and mucking. Support is not considered, because of the good quality of the rock. The rock mass rating (RMR) geomechanical classification of the rock mass [13] carried out by MSP indicates that the rock mass in the tunnel section can be classified as very good rock, with RMR over 85 points (class 1), and so no support is required.
After measurements had been made of the operational times of the cycle in the exploratory activities with characteristics similar to those that will be simulated, under the guidance of the Engineering Department of MSP, each of the activities was characterized in a statistical analysis that allowed determination of the probability distribution that best fitted the performance of each of the unit operations. Table 2 shows a summary of the statistical analysis, with the corresponding assigned PDFs.

Ventilation was kept constant in time, because MSP provides lunch for workers at the same time, and the duration of both lunch and ventilation is 90 min.
The distributions shown in Table 2 are those that were used to produce the algorithm and will be used in the Monte Carlo simulation of the tunnel construction time.
In general, the analysis presented in Table 2 is based on data obtained from field sampling, and the fitting was made by MSP, who are solely responsible for the data handling, but it should be noted that the fitting of the probability distribution curves was done by the AndersonDarling method.
The simulation is involved between 10^{5} and 10^{6} iterations. It was found that the variability between the first and last simulations from this interval was not significant, considering the mean and the mode of the results (see Table 3), so it is believed that all the results delivered by simulation beyond 10^{5} iterations are good.

As stated in the model, every simulated event involves the construction of a tunnel with the specified length, taking into account the duration of the shift and the tolerance (“Tolerance” in Table 3) to ending the activities within the cycle. This means that less time is left rather than stipulated in this item to end the shift, so the activities are suspended and are taken up by the following shift.
Finally, it is also necessary to consider the restarting time (“Reset” in Table 3) of the activities, which corresponds to the time required to resume activities if the operation is interrupted by the end of the shift.
Both the tolerance and the restarting time can be modeled as Markov chain processes [1], because the success or failure of an event has an effect on the following one.
For the simulation, a tunnel of 3.5 m × 3.0 m with a length of 560 m is considered, with two work shifts per day, each of 540 min duration, with a tolerance of 60 min, which represents the time between the end of a cycle and the end of the shift. If the tolerance time is less than 60 min, the activities are finished and service or other activities related to the operation, such as cleaning, are carried out.
A 30 min time is considered for resumption of activities, implying that if the cycle time is longer than the shift time, this value is added to the cycle time, because all the distributions presented in Table 1 consider the starting time of the activity, but not a restart caused by a shift change. The considered values correspond to expert data coming from the experience of MSP in the operations area.
8. Analysis of the Results
According to the data presented in Table 3, and considering a simulation with 10^{5} iterations (the difference compared with a simulation with 10^{6} iterations is negligible and the shorter simulation is easier to handle with the available statistical software), Figure 2 shows the probability distribution obtained for the simulation in Mina Romero.
In contrast to the conventional planning method, which determines one value for the required number of shifts, one of the advantages offered by the simulation method presented here is that it is possible to have both pessimistic and optimistic scenarios with respect to the number of shifts required for carrying out the work. These scenarios can be considered as the lower and upper limits of the confidence interval that describes the construction time for the tunnel under study in shifts.
The simulation produces a histogram with a mean of 266.89 shifts, a mode of 267 shifts, a median of 267 shifts (with a minimum of 257 and a maximum of 277), and a standard deviation of 2.31, with a distribution that is symmetric in form. This histogram is fitted into a gamma distribution (Figure 2) with parameters , , (minimum), and (maximum). The curve was fitted using the AndersonDarling method, as mentioned earlier.
For the case in question, 267 shifts, which is the value of the mean as well as the median and the mode, were considered. It was decided to use this value because it is a good representation of the simulated case and is the most repeated value. Figure 3 shows the probability distribution for the simulation in Mina Romero, and from this curve, it is possible to obtain the probability of success, in this case 0.6 (60%).
The gamma distribution is not of symmetrical type, but in this case is the best fit to the simulation data using the AndersonDarling method—this is why the mean probability of success is 60% rather than 50%. Figure 2 shows that the mean is not in the middle of the curve, and this is confirmed by Figure 3.
Once the tunnel construction, which took 133 days, was finished, the means resulting from the simulation could be compared with the MSP plans and with the actual data (Table 4), revealing an error of 0.37% compared with the real time, while with the conventional planning method used by MSP the error was 12.28% with respect to the actual construction time. This difference is significant because, when it is translated into execution days, it can be seen that the conventional planning method used by MSP underestimated by 16.5 days the time required, whereas the simulation gave a mean differing by only one day from the real execution time. When the standard deviation is considered, we have a quite precise tool for planning, because the value considered is 267 ± 3 shifts.

It should be mentioned that the plan made by the Engineering Department of MSP is not represented in the histogram shown here. A more careful analysis would show that there is no event similar to what was planned and indeed that it would be very difficult for the event planned by the mining company to occur.
As can be seen, this simulation methodology based on the Monte Carlo method can estimate the time required for the construction of a tunnel used in underground mining.
9. Conclusions
The data analysis performed at MSP shows that this kind of stochastic simulation is a very effective tool for planning the construction time of a tunnel.
Beyond the accuracy of the means, the range of minima and maxima obtained by the simulation is interesting, because it delivers a potentially useful parameter for establishing planning criteria.
Based on the minimum and maximum values obtained from the simulation, optimistic and/or pessimistic scenarios can be proposed that they can serve as background information for making mine planning decisions.
Because of the random nature of the execution times of the operations involved in mining construction, it has been determined that a planning methodology based on the Monte Carlo method provides a better fit to real conditions than a conventional approach, because, with its use of probability distributions, it incorporates the variability inherent in the planning process.
By incorporating variability into planning, it is possible to determine with greater certainty the ranges over which the execution times of the different operations fluctuate. It is thereby possible to reduce the financial risks due to planning errors while maximizing the exploitation of resources.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported by DICYT/Universidad de Santiago de Chile, USACH (code 051515VN), and PPGE3M Universidade Federal do Rio Grande do Sul.
References
 F. S. Hillier and G. Lieberman, “Markov chain,” in Introduction to Operations Research, pp. 675–682, McGrawHill, 9th edition, 2010. View at: Google Scholar
 S. P. Singh and P. Xavier, “Causes, impact and control of overbreak in underground excavations,” Tunnelling and Underground Space Technology, vol. 20, no. 1, pp. 63–71, 2005. View at: Publisher Site  Google Scholar
 F. T. Suorineni, P. K. Kaiser, and J. G. Henning, “Safe rapid drifting—support selection,” Tunnelling and Underground Space Technology, vol. 23, no. 6, pp. 682–699, 2008. View at: Publisher Site  Google Scholar
 Gobierno de Chile, Código del Trabajo, Capítulo IV ‘De la jornada de trabajo’, Gobierno de Chile, 2012.
 N. Metropolis and S. Ulam, “The Monte Carlo method,” Journal of the American Statistical Association, vol. 44, no. 247, pp. 335–341, 1949. View at: Publisher Site  Google Scholar
 H. T. Chiwaye and T. R. Stacey, “A comparison of limit equilibrium and numerical modelling approaches to risk analysis for open pit mining,” Journal of the Southern African Institute of Mining and Metallurgy, vol. 110, no. 10, pp. 571–580, 2010. View at: Google Scholar
 E. Ghasemi, K. Shahriar, M. Sharifzadeh, and H. Hashemolhosseini, “Quantifying the uncertainty of pillar safety factor by Monte Carlo simulation—a case study,” Archives of Mining Sciences, vol. 55, no. 3, pp. 623–635, 2010. View at: Google Scholar
 R. Khalokakaie, P. A. Dowd, and R. J. Fowell, “Incorporation of slope design into optimal pit design algorithms,” Transactions of the Institution of Mining and Metallurgy Section A: Mining Technology, vol. 109, no. 2, pp. A70–A76, 2000. View at: Publisher Site  Google Scholar
 M. A. Morin and F. Ficarazzo, “Monte Carlo simulation as a tool to predict blasting fragmentation based on the KuzRam model,” Computers & Geosciences, vol. 32, no. 3, pp. 352–359, 2006. View at: Publisher Site  Google Scholar
 M. Sari, C. Karpuz, and C. Ayday, “Estimating rock mass properties using Monte Carlo simulation: ankara andesites,” Computers and Geosciences, vol. 36, no. 7, pp. 959–969, 2010. View at: Publisher Site  Google Scholar
 I. M. Sobol, A Primer for the Monte Carlo Method, CRC Press, Boca Raton, Fla, USA, 1994. View at: MathSciNet
 X. Chen, “Resuing shrinkage stoping,” Engineering & Mining Journal, vol. 199, no. 10, pp. 34–36, 1998. View at: Google Scholar
 Z. R. Bieniawski, “Rock mass classification in rock engineering,” in Proceedings of the Exploration for Rock Engineering, pp. 97–106, A.A. Balkema, Cape Town, South Africa, 1976. View at: Google Scholar
Copyright
Copyright © 2015 Juan P. Vargas et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.