Discontinuum Modelling in Rock Engineering

. The evolution of numerical modelling in rock mechanics has made possible a more realistic representation of the behaviour of rock masses, and a more reliable prediction of the response of rock engineering structures. Discontinuum modelling techniques, namely the discrete element method, based on the explicit representation of the rock mass discontinuous structure, have progressively acquired a broader role. In this lecture, the essential concepts and main features of these numerical techniques are discussed, with reference to their historical development. The use of these idealizations is exemplified by two areas in which they proved very effective. On one hand, the analysis of practical rock engineering problems, intended, for example, for safety assessment or monitoring interpretation. This type of application is illustrated herein by dam foundation analyses based on deformable block models. A second level of discontinuum modelling involves research aimed at the understanding of the fundamental behaviour of rock and rock masses. An example of this approach is the analysis of fracture phenomena with particle models, primarily of laboratory specimens, but being extended to larger scales. Critical issues identified in the effective application of these modelling tools are examined, as well as some foreseeable evolution trends.


Introduction
The development of Rock Mechanics as an independent discipline, with its own field of knowledge and methodologies, may be associated, to a large extent, with the progressive recognition of the role played by discontinuities in the behaviour rock masses. In 1964, in a well-known paper presented to the 8th ICOLD Congress in Edinburgh, Manuel Rocha stated that among the difficulties then facing Rock Mechanics should be mentioned "in first place, the consideration of families of joints and faults which turn rock masses into systems of more or less adjusted blocks" (Rocha, 1964a). Then, he questioned "how far such systems can be treated by extending to them the laws of behaviour applying to media assumed to be continuous or if it is necessary, at least in certain cases, to develop a mechanics of media divided by families of discontinuities". The examination of elementary systems, such as those depicted in Fig. 1, was a step towards the definition of appropriate conceptual models. In the same paper, Manuel Rocha recognized the limitations of the available analytical methods to deal with systems of such nature and complexity, pointing to physical experiments as a promising means to simulate discontinuous media and provide understanding about their mechanics.
In the late 1960s, numerical methods were advancing into the representation of systems cut by discontinuity surfaces, whether within the framework of the finite element methods, or in the newly developed discrete element method, a path that will be discussed in more detail in the following section. However, it would take a decade for the digital computing hardware to be able to handle large, complex 2D models, and perhaps two decades for 3D representations of problems of practical interest to be feasible. Today, we have the means to analyse problems with large sizes and intricate joint patterns. However, the hardware limits continue to be challenged, for example in the study of dynamic processes or coupled physics problems. In fact, the range of physical phenomena that modelling is called upon to simulate has significantly enlarged, with new requests arising in areas such as waste isolation, oil and gas exploration or CO 2 sequestration. Fluid flow has always been recognized as a major factor in rock engineering, namely in dam foundation problems, also being a problem in which discontinuities, and the way in which they are connected, play a dominant role. Transport and thermal phenomena have also to be addressed in some of the expanding fields. The need to consider the coupling between all of these processes adds to the difficulty in model development, and, more critically, in their validation for reliable application.
role they assumed in the study of fracture phenomena. The need to examine fracture processes in rock based on a consistent analytical framework, while taking into consideration its specific aspects arising from heterogeneity and grain structure, has long been felt in rock mechanics (Fairhurst, 1972). The development of discrete element models based on circular or spherical particles has made numerical models a key ingredient in this important field. The power revealed in the analysis of laboratory test specimens is prompting their use at larger scales, which is again pushing the limits of computer technology.
The present lecture is aimed at the discussion of the essential issues that presently face the application of discontinuum models in rock engineering. The underlying conceptual models will be first examined. Then, it is instructive to look into the historical development and evolution of the numerical models, in particular those based on the discrete element method. Being impossible to cover all areas of rock mechanics, two fields of application were selected to discuss the use and capabilities of these models, one related to engineering practice, the other more research oriented. First, the problem of dam foundations in rock will be addressed, an application of block models at the engineering scale, which involves specific aspects such as hydromechanical analysis, seismic analysis and failure assessment methods. Secondly, the analysis of fracture processes by means of particle models is examined, in particular looking into the potential of these approaches for broader problems, as they begin to extend from the research environment to the engineering office. Finally, some concluding remarks are presented on the current trends and development needs of discontinuum modelling in rock engineering.

Representation of a rock mass from an engineering perspective
An engineering model is necessarily a simplification of the physical reality, more often intended to answer a spe-cific question, e.g. about safety or performance of a proposed design or an existing structure, rather than to provide a meticulous description of nature. The influence of the rock discontinuities on the overall response of the rock mass depends on the specific questions and behaviour traits being investigated. Therefore, a given geological condition may be idealized in different ways depending on the problem at hand. Whether the purpose of the analysis is preliminary or detailed design, whether the aim is the interpretation of monitoring data or the safety reassessment of an existing structure may have a strong influence on the choice of a conceptual model and its numerical implementation, as different geological or geomechanical traits assume more or less significance. It is also likely that the knowledge about the rock mass characteristics differs in each case, as well as the level of accuracy sought in the solution, or the time available to obtain the results. All of these may call for different ways to idealize the problem, giving more weight to the representation of the features recognized as dominant.
Two fundamental options exist for the representation of a jointed or fractured rock mass. The essential difference lies in the manner of representing the discontinuituies (faults, joints, fractures, interfaces) and their effects on rock mass behaviour. On one hand, we have the equivalent continuum approach, in which a continuum constitutive model is employed to represent in an average manner the effects of the discontinuities. The alternative is the discontinuum approach, in which the discontinuity surfaces are explicitly represented individually. The rock mass is then viewed as a jointed medium, or as an assembly of interacting blocks (or particles).
One of the basic criteria in deciding whether to use an equivalent continuum or blocky model is based on the scale of the jointing with respect to the size of the excavation or foundation under study, as illustrated by the classical scheme by Wittke (1990), reproduced in Fig. 2. Closely spaced joint sets are more appropriately included in the rock matrix behaviour, with explicit representation employed only for major, singular features. Major faults at known locations typically have to be considered individually. In practice, a range of representation options are possible combining continuum and discontinuum in the way that best fits the needs of the problem at hand. For example, in failure analysis, it may be sufficient to include a few joint planes of each set to represent the major potential mechanisms. However, the deformability of the block material, if it is important, should include the contribution of the joints not explicitly modelled.
The equivalent continuum option simplifies the geometric definition, but requires more complex material models. In some cases, when the loading is not expected to cause non-reversible effects, elastic models may be used, possibly with anisotropy to account for the joint orientation. Early no tension material models were  (Rocha, 1964a). followed by elasto-plastic and damage models (e.g., for conditions of very fractured rock with statistical isotropy), and then by multilaminate or smeared crack models, capable of providing a directional variation of tensile and shear strength. In problems where the deformability of the rock mass is the main issue, or the profusion and complexity of joint patterns may be difficult to simulate by a block pattern, these continuum representations are a good option. The well-known discussion of the range of application of the Hoek-Brown failure criterion illustrates the continuum vs. discontinuum conceptual divide (Hoek & Brown, 1980).
The discontinuum option requires more geometric information, more detail, but it generally involves simpler constitutive laws. Jointed or blocky systems tend to represent more faithfully the types of experimental response observed in the field. However, given the natural uncertainty about jointing patterns, it may always be questioned whether a selected joint structure is representative. A relatively crude representation of the real jointing patterns in a discontinuum model may be justified if it addresses the essential impact it has on the response. The simple idealizations discussed by Rocha (1964a), shown in Fig. 1, illustrate key aspects that govern behaviour, such as orientation and persistence of each joint set or imbrication of crossjoints. Current research on discrete fracture networks is providing new tools to generate complex patterns, whenever there is statistical data to support them.

Discontinuum numerical models. The discrete element method
In the 1960s, the newly developed finite element (FE) method became an important tool in rock mechanics, as it allowed the stress analysis of complex geometries which were not addressable by the existing closed form solutions. FE models provided the means for the equivalent continuum representation of rock masses, initially based on linear elasticity, but soon afterwards extended to more elaborate non-elastic constitutive relations. The need of a truly discontinuous numerical representation, however, remained unanswered, particularly when failure mechanisms in jointed or blocky rock had to be assessed. The first solution to this problem was the joint element proposed by Goodman et al. (1968), which allowed the explicit representation of individual discontinuities in a FE model of a rock mass. This is a special type of finite element, which is assumed to have zero thickness and employs a joint constitutive model relating the joint normal and shear stresses with the differences in nodal displacements across the element, unlike the stress-strain laws of continuum elements. In the elastic range, therefore, the joint normal and shear stiffness parameters, as defined in rock mechanics, determine the element deformability. Non-linear behaviour in tension or shear was included from the outset, while extension to 3D posed no theoretical difficulties, only demanded more computer power. Finite element models with joint elements were capable of actually representing a discontinuous medium, realizing the conceptual model of a continuum crossed by a few discontinuities where failure may take place. There were, however, some limitations when dealing with complex jointing patterns involving multiple sets leading to blocky systems. For example, failure mechanisms involving block rotational modes, with large movements, were not accounted for, and the small displacement assumption, prevalent in most codes, prevented changes in the system connectivity. In addition, the numerical methods of solution in use at the time, based on iterative methods using stiffness matrices, did not permit total block separation. These limitations were particularly obvious in slope stability analysis, and it was to address these problems that Cundall proposed a new numerical approach, designated as the 'Distinct (or discrete) element method' (DEM) (Cundall, 1971). This method viewed the rock mass as an assembly of distinct blocks or particles in mechanical interaction. These blocks could behave as rigid bodies, as often assumed by conceptual models of rock slope stability, a simplification that was not possible with the existing finite element models. Toppling failure was one of the first applications that illustrated the capabilities of the new method (Fig. 3).
Earlier concepts of blocky systems, such as the "clastic mechanics" proposal of Trollope (1968), had a limited practical application given the shortcomings of the analytical framework. The key to the success of Cundall's DEM approach was the solution method chosen, which relied on the integration of the equations of the motion of the blocks using a time stepping explicit algorithm (Fig. 4). This type of solution method was well known in the finite difference community, as dynamic relaxation was one of the techniques employed to solve systems of linear equations before computers made direct solution feasible (e.g., Otter, 1965). Physically it amounts to seeking the static solution of a mechanical system by means of artificially high damping. It is particularly suited for discrete systems undergoing large motions, for example, in the analysis of failure mechanisms, disaggregation processes or rock falls. Programming the contact detection and update algorithm in an efficient manner became a crucial component of discrete element codes, freeing the user from the need to define manually the initial system connectivity and the way it evolves during the analysis (e.g., Cundall, 1988).
A distinctive feature of discrete element models is the representation of contact, based on sets of contact points. Rock joints are treated numerically as surfaces of interaction between neighbouring bodies, instead of being predefined as joint elements. This greatly enhances the flexibility of model generation and discretization, and also simplifies the consideration of large movements (e.g., Lemos, 2008). Assigning representative contact areas to these contacts allows the standard rock joint constitutive models to be employed, such as the classical Mohr-Coulomb model, or more elaborate ones, such as Barton's joint model (Barton et al., 1985).
In order to be applicable to wider problems, block deformability had to be considered. First, the assumption of uniform stress field in each block was adopted, but its limitations were evident, as complex deformation modes, such as bending, could not be represented. The alternative was to divide each block in an internal finite element mesh, which could be refined according to the stress analysis precision required. This approach to "deformable blocks", developed  in 1978, was used in the new UDEC code (Cundall, 1980), and its 3D counterpart 3DEC (Hart et al., 1988). Similar schemes, sometimes designated as "discrete-finite elements", are now employed in various codes. In Cundall's codes, the solution method for deformable block models remains explicit, with the degrees of freedom of the system now being the displacement of the nodal points, instead of the rigid block displacements and rotations.
In the original 1971 paper, Cundall already included the circular particle models, treated numerically in the same way as rigid blocks. However, the geometric calculations involved in detecting and updating contact between circular (or spherical) particles are much simpler than between polygonal (or polyhedral) blocks. This implies that much larger systems can be simulated with reasonable run times, a major factor in the fast expansion of particle models. Initially proposed for micro-mechanical modelling of soils (Cundall & Strack, 1979), they became a choice tool for rock fracture studies.
Presently, many DEM codes, based either on block or particles, are routinely employed in rock mechanics research and practice, as well as other numerical methods that share some of its essential concepts for the representation of a discontinuum, including Discontinuous Deformation Analysis (DDA) (Shi & Goodman, 1988), Discrete-Finite Element Method (DFEM) (Munjiza, 2004), Non-Smooth Contact Dyanmics (NSCD) and others (e.g. Jing & Stephansson, 2007;Zhao et al., 2012).

The use of discrete element models and codes
The ability to analyse systems of discrete bodies, in particular after the development of efficient techniques to cope with very large systems, has turned the discrete element method into a powerful instrument for detailed modelling in many fields of physics and engineering. Closer to rock mechanics are various applications involving geomaterials, such as soil, rockfill, masonry or concrete. Analysis of rock falls, mine caving and fragmentation processes are also problems for which DEM is particularly suited. A particularly intensive research effort has been directed towards addressing broader physical problems, for example involving thermo-hydro-mechanical and transport processes. This broad field of application has led to a multitude of formulations and software.
The variety and sophistication of the numerical techniques available, as well as the codes in which they are implemented, often present difficulties to the student or practitioner in the choice of the most effective tool for the problem at hand. A clear understanding of the underlying assumptions of a given code, in order to judge its capability to reproduce the conceptual mode adopted, is the first requirement for proficient use. The development of proper modelling methodologies of application of numerical models needs to receive more attention, particularly if we consider the natural tendency towards more complex models.
The amount of detail that may be included in a model is always limited by the experimental information available. Often, the choice of a continuum model is simply dictated by the preliminary nature of the study. Starfield & Cundall (1988) presented an influential discussion of analysis methodologies in rock mechanics which stresses the role of models in providing understanding about a given rock engineering problem. A numerical model may be viewed as a numerical laboratory in which to test different hypotheses of behaviour by comparison with the available knowledge, and to identify the critical factors, e.g., those parameters that need reliable experimental evidence since they play a dominant role. Analysis of a given problem may proceed by a sequence of models intended to clarify specific aspects of behaviour. To address the entire range of issues with a single very detailed representation is often inefficient. Calibration of a multitude of parameters in a complex model based on a few known response features can be very misleading if the physical role of each one is not well understood.
The ability to automate model generation or output processing tasks increases the power of the simulation tools, but demands more expertise. Robust methodologies for model building and verification are essential for effective engineering analysis. Facilities to implement userdefined rock or joint constitutive models are being offered by advanced codes. In fact, versatility is a highly desirable quality in engineering software, allowing users to progress from elementary or occasional use to elaborate models, always building on the acquired expertise. In addition, it is of crucial importance to be proficient in the representation, processing and interpretation of the numerical results, so that as much knowledge as possible may be gained from the analysis.

Analysis needs and modelling choices
The design of concrete dams is subject to the fulfilment of general requirements of serviceability and safety shared by all engineering structures. The former imply that the structure should be able to perform its function satisfactorily, under the normal operating conditions. The latter require that the structure should not create unacceptable risks for those who use it or the public. In the case of dams, it is mandatory that, even in extreme conditions, uncontrolled loss of the retained water does not occur (Pedro, 1995).
The verification of the two types of requirements entails different analysis needs. For normal operating conditions, the aim is that the behaviour of structure and foundation remain essentially in the elastic, reversible domain. Therefore, equivalent continuum models are usually sufficient to represent the rock mass. The main requirement is to represent correctly the deformability of the various rock mass zones, and possibly major singular features. These models are typically sufficient to predict structural displacements during normal operation, and may be validated and calibrated with monitoring data.
Safety assessment involves the analysis of the possible failure modes of the structure or foundation, considering extreme load conditions, including flood water levels and large earthquakes. Failure mechanisms involving the dam-rock interface and rock mass discontinuities need to be considered, for which discontinuum conceptual models are particularly appropriate. Rocha (1978) discusses the possible failure modes of gravity dams (Fig. 5), which in the simplest cases may be resolved by analytical calculations. Numerical models are now the best tools for these analyses, having no difficulty in addressing, for example, the case of failure surfaces with multiple segments.
The accident of Malpasset brought into attention the role of the hydraulic behaviour of the rock mass, and the importance of water pressures on rock discontinuities. In the 1960s, analytical and graphical methods were developed to assess the stability of rock wedges in arch dam abutments, including uplift forces, as shown in Fig. 6 (Londe, 1973).
Physical models were for a long time the main tool for stress analysis of arch dams, which only permitted the deformability of the foundation to be represented, and in a rather simplified manner (Rocha, 1964b). Experiments with blocky foundations were also carried out to simulate more complex rock mass conditions (Oberti & Fumagalli, 1963). Sliding failure mechanisms may also be analysed, as in the case of an arch dam studied at LNEC and depicted in Fig. 7 (Gomes, 2006). Numerical models gradually became the main analysis tool, but the contribution of physical models to their validation should not be overlooked. For the case of dynamic slip, shaking table tests remain an important source of data, given the complexity and variability of these phenomena.
Earlier numerical models were 2D representations or gravity dams adopting a continuous foundation. In spite of being limited by linear elastic assumptions, the numerical models had the means to include the spatial variability of deformability and anisotropy, a significant step forward from a practical point of view. Discontinuous finite element models using Goodman's joint model were soon applied to simulate major discontinuities. The discussion by Rocha (1978) of the analysis of Água Vermelha dam, in Brasil (Fig. 8), in which non-linear joint elements were used to represent a horizontal discontinuity in the foundation provides a clear example of the insight into the system behaviour that could be gained from a pertinent use of the new numerical techniques. 142 Soils and Rocks, São Paulo, 36 (2)   The failure mechanisms illustrated in Figs. 4 and 5 may now be straightforwardly approached by DEM block models, which makes them a very helpful tool in dam design. As seepage and uplift pressures are of critical importance on performance and safety issues, they must be properly taken into account in the analyses.

Hydromechanical behaviour
The significance of uplift pressures for dam safety was recognized in the pioneering work of Lévy, following the accident of Bouzey dam in 1895, which made clear some inadequacies in traditional designs and the importance of drainage (e.g., Bretas et al., 2011). The failure of Malpasset Dam in 1959 at the completion of the first filling of the reservoir, as a result of sliding on a rock wedge under the left abutment, originated a great research effort aimed at a better understanding of the behaviour of rock foundations, in order to improve design practices (Duffaut, 2011). Londe's analysis stresses the influence of the jointing and deformation patterns on the rock mass permeability, which allowed high water pressures to be installed and favoured the sliding mechanism (Londe, 1987).
The analysis of the water flow and head distributions under dams was a challenge for a long time. Electrical analogs and other physical experiments were tried, before numerical solutions became the standard to study flow governed by Darcy's law. Finite-difference methods were initially used (e.g., Serafim, 1968), followed by finite element models, first assuming uncoupled flow, and later coupled flow-stress analysis. Cases of gravity and arch dams are presented, for instance, in Wittke (1990). Multilaminate models were also used to account for the joint structure (Lamas & Sousa, 1993). The experimental study of flow in fractures by various researchers provided the background for discontinuum numerical representations, either using joint finite elements, or within a DEM framework (Damjanac & Fairhurst, 2000). Several authors have applied 2D fracture flow models to the study of the water flow under gravity (e.g., Lemos, 1999;Gimenes & Fernández, 2006). Whether joint patterns are idealised, or some type of random generation is employed (Barla et al., 2004;Fig. 9), the explicit inclusion of joints planes is particularly important for stability assessments.
While there is no question that flow in a rock mass is more realistically represented by fracture flow models, a number of difficulties often hinder the practical viability of this approach. First of all, fracture flow models require much more information than a continuum idealization. For the latter, the average permeability of the various foundation zones is sufficient, while a discontinuous representation of flow implies the specification of joint apertures, as well as realistic values of the joint normal stiffness since this parameter governs the stress-flow coupling. These data may not be available. The representation of the joint network also requires a good knowledge about the rock mass structure. For dams, the critical areas for flow are often at Soils and Rocks, São Paulo, 36 (2) (Rocha, 1978). shallow depths below the dam, typically more disturbed and with higher variability. Furthermore, the representation of the grout curtain is also more problematic with a fracture flow model. A study of the hydraulic behaviour of Alqueva dam, involving packer tests to detect locations of water inflow into drainage boreholes, has clearly shown the importance of local conditions on the flow (Farinha et al., 2011). A continuum flow model (Fig. 10), where the grout curtain and drainage systems were represented, was calibrated with the in situ tests data, but a detailed fracture flow model would not be viable, particularly at the design stage, when monitoring results were not yet available.
As discussed in the next section, discontinuum models are the ideal tool for failure analysis, given the dependence of the most likely collapse mechanisms on the rock mass features. In these safety assessment calculations, water pressure fields need to be specified in all discontinuities, so that the possibility of sliding is properly checked on the basis of effective stresses. However, these water pressures need not be the result of a fracture flow analysis, if there is not enough data to perform one. Either standard design hypotheses are assumed, or a less demanding continuum flow calculation may be performed, from which the water pressures can be transferred to the mechanical model for safety assessment. It should also be stressed that a fracture flow model requires a much finer representation of the jointing, in particular its connectivity, while for a stability analysis only the main discontinuities that may contribute to define a failure mechanism are necessary.

Failure assessment
Analysis of failure of gravity dams has been traditionally performed by 2D models, as upstream-downstream sliding is the major concern. The dam-rock interface is the first location to be checked, followed by other mechanisms defined by the rock joints (Fig. 5). Sliding mechanisms involving multiple segments exceed the capabilities of the classical analytical calculations, so numerical models become the best option. Rock discontinuities have obviously a 3D structure, which may have to be considered in stability analysis, namely for blocks supported in the valley slopes (e.g., Lombardi, 2007). The water pressure distribution is often the key parameter, stressing the importance of effective drainage. The strength properties of the concrete-rock interface encompass some uncertainty, so many code regulations require the consideration of a friction-only hypothesis (e.g., ICOLD, 2004).
For arch dams, Londe's classical analysis (Fig. 6) illustrates the conceptual model of abutment failure. Numerical representation may be performed with finite element models with joint elements (e.g., Wittke, 1990;Alonso et al., 1994). At present, DEM codes provide the most flexible option to conduct this type of evaluation. Fig. 11 shows a model of Baixo Sabor dam, a 130 m high concrete arch in a wide valley in the North of Portugal (Matos & Paixão, 2007). The model, created with 3DEC (Itasca, 2007), may be used to exemplify the main issues posed by this type of analysis (Lemos, 2012b). All the blocks, in the dam and foundation, are deformable and assumed elastic. The dam is divided into cantilevers defined by the contraction joints, each one represented by means of higher order finite elements, with an accurate performance in bending. For the rock mass, deformable polyhedral blocks are employed, discretized in tetrahedral element internal meshes, which allow a proper consideration of the spatial variation of the rock mass deformability.
In this model, a very simplified representation of the rock mass discontinuities was chosen, focusing on the area of concern, the right abutment. First, the major faults were included at their known locations. Then, the block system was created by considering only a few joint planes of each of the 3 main sets of the granitic rock mass, in such a way that the various foreseeable failure modes could be generated. Away from this area of interest the model was ex- tended up to the boundaries by means of large elastic blocks. In the absence of detailed information to perform a flow analsys, conservative assumptions were considered to assign water pressure along the dam-rock interface and in all the rock joints. The safety quantification was performed by reducing the frictional strength in all discontinuities until a failure mechanism developed, as depicted in Fig. 12. The evolution of various displacement indicators throughout the strength reduction process, shown in Fig. 13, confirms the ample safety margin, as significant movements only take place for reduction factors above 2.
The question of how much complexity may be added to a workable model is certainly subjective, and different approaches are feasible. Some engineers prefer to use simpler models to study separately specific zones of concern. Different representations of the same dam may be em-ployed, for example to study right and left bank stability, or different joint set locations or combinations. Others prefer to build a single large model in which all possible failure modes can be generated. This implies extra time in model preparation and verification, and the interpretation of the results may be less clear or straightforward. In fact, models have the power to provide insight into the traits of structural behaviour, which is sometimes more useful than yielding a safety factor number.
In the model shown above, the rock mass discontinuities were represented by extended planes, a justified assumption for the existing faults, but clearly conservative for the other rock joints needed to create a failure mode. In this case, a comfortable safety factor was achieved, even neglecting the effects of rock bridges. For other cases, this assumption may be too conservative, so a more elaborate Soils and Rocks, São Paulo, 36 (2)  representation of the jointing would be desirable. Generation of discrete fracture networks for use in stability analysis remains an important research topic, in order to improve the realism of rock mass representations, in agreement with its statistical descriptions derived from characterization studies (e.g., Lorig, 2007). However, for analyses with non-persistent joints to be dependable, they need to take into account the possible breakage of the rock bridges and the progression of the failure surface. This issue will certainly gain from research on detailed fracture models, to be discussed below.

Seismic analysis
Seismic analysis of concrete dams traditionally concentrates on the response of the structure to assess the possible occurrence of cracking or collapse, with the rock mass being represented as an elastic continuum. For gravity dams, analysis of sliding on the foundation joint under seismic loading is now often performed for large design earthquakes, in order to provide estimates of permanent displacements. Oversimplified models in which the dam is represented as a single rigid block are not advisable, since the dam dynamic deformation alters the stress distribution on the foundation plane and thus the slip movements. Postseismic stability analysis, considering possible damage to the dam-rock interface and to the grout curtain, with its potential implications in increased uplift diagrams, is also an important use of numerical models (e.g., Aillard & Léger, 2008).
When non-linear behaviour of the rock foundation needs to be considered, more elaborate numerical setups are required, as the classical massless foundation hypothesis is no longer acceptable. The inertial and dissipation properties of the rock mass have to be included, with the appropriate non-reflecting and free-field boundary conditions (e.g., Lemos, 1999). Seismic analysis of arch dam foundations is a topic that has been seldom addressed in the past, but the tools are now available to study the dynamic response of complex block systems. An example, using the model of the previous section, is shown in Fig. 14 (Lemos, 2012b). The permanent displacements after an earthquake with a peak ground acceleration of 0.5 g, depicted in the figure, are not expected to pose any risk to the dam integrity.
The transient dynamic water pressures that can be generated in a rock joint have been investigated experimentally (Javanmardi et al., 2005), but their representation in numerical models still poses some difficulties. The simpler option is to disregard the transient pressure effects, keeping the hydrostatic water pressures unchanged, as it was done in the model of Fig. 13, which may be a conservative simplification in many cases (e.g., Lemos, 1999). Given the uncertainties about these phenomena and the dispersion of estimates of dynamic slip that are usually obtained with numerical models, as well as shaking table tests, multiple analyses are always indispensable, varying the seismic input records and the less well-known properties.

Detailed DEM models
The discontinuous nature of geo-materials is evident at various scales, for each one being possible to identify different dominant phenomena and behaviour patterns. At the scale of engineering structures, when global stability is the main concern, the large block systems analysed in the previous section address the representation of the rock mass major features. Closer inspection would reveal a network of fissures which govern, for example, the patterns of fluid flow. Below that scale, the rock matrix may be described by the grain structure and its micro-cracks, either inter-or intra-granular. This is sometimes referred to as the mesoscale. Further down, we might consider the molecular scale, but that is beyond the scope of rock mechanics. While the definition of the entities and the modes in which they interact will vary, discrete elements are sufficiently versatile for the numerical treatment of all of these problems. Polygonal blocks are a natural representation of jointed rock. At the finer scales, a higher degree of idealization may be adopted, in terms of simpler particles shapes. Circular particles are especially efficient numerically, since the geometric resolution of the contact is straightforward, therefore allowing very large systems to be assembled. The representation of the behaviour of complex systems starting from elementary interactions is common to other scientific fields. In physics, the techniques of Molecular Dynamics (e.g., Pöschel & Schwager, 2005) are based on concepts similar to the DEM. As computer power increases, the ability of these methods to address problems with practical interest has grown significantly, and their application has expanded. Cundall's original code (1971), already had circular particles in addition to blocks (Fig. 3), but the complete formulation of circular particle models was established by Cundall & Strack (1979). The leading motivation at the time was the micro-mechanical study of soils based on a conceptual representation of granular media as assemblies of rigid circular particles. DEM was employed to investigate the fundamental mechanical response of randomly generated assemblies of particles, as a tool to understand the underlying elementary phenomena, in order to allow the investigation and testing of constitutive assumptions, which could be then employed in continuum constitutive laws.
However, it is the study of rock fracture that has become the key application of circular particle models in rock mechanics. By connecting the particles with bonds, which are allowed to break according to specific failure criteria, the random assemblies of disks or spheres provide a powerful representation of the rock matrix, simulating the complexity of the fracture propagation paths. While circular particles have the advantage of simplicity and speed, polygonal blocks provide a more realistic representation of the interlocking grain structure, so they are increasingly resorted to in research work.
The next 3 sections will focus on the application of bonded particle models to represent, respectively, the rock matrix behaviour, rock joints, and, finally, the comprehensive response of rock mass.

Essential concepts
The well-known paper by Potyondy & Cundall (2004) marked a decisive step in the numerical study of rock fracture by DEM, and is an excellent introduction to the designated Bonded Particle Models (BPM). The basic idea is to represent the rock material as an assembly of disks (in 2D) or spheres (in 3D), randomly generated according to some size distribution curve. These rigid particles are linked by breakable bonds. The system deformation derives from the inter-particle movements. Very few property parameters are necessary. The contact normal and shear stiffness govern the deformability of the unbroken system. Bond strength is described by the tensile strength and cohesion terms. After bond breakage the interaction between the particles becomes purely frictional. The bond formulation employed, known as parallel bonds, offers the capability to transmit moments, in addition to forces, between the particles.
This very simple set of assumptions allows very complex forms of behaviour to arise as the random particle assemblies are progressively loaded, which can be compared with experimental evidence. Fig. 15 shows results of 2 biaxial tests, the first with a low confining stress of 0.1 MPa, and the second with a high confining stress of 70 MPa. The expected axial splitting mode of failure is obtained in unconfined tests, while shear bands develop in the confined tests, with more distributed damage (represented by breakage of bonds) as the confining stresses increases. The transition from brittle to more ductile behaviour is thus simulated by the numerical model. However, this model does not reproduce correctly some aspects of experimental behaviour, namely the ratio of uniaxial tensile to compressive strength is too high, and the macroscopic friction angle, inferred from biaxial tests at different confining stresses, is too low. Clearly, the friction interlocking present in a real rock matrix is not being simulated. The various procedures since developed by various researchers to overcome these deficiencies will be discussed in the following section. Several important aspects about this approach should be noted. First, the importance of the procedure to create the random assembly, discussed in detail by Potyondy & Cundall (2004), which is required to provide uniform, isotropic samples with the intended porosity and size distribution. Secondly, a major issue involved in this approach is the determination of the micro-properties, namely the bond stiffness and strength parameters, which in general cannot be directly related to the desired macroscopic properties. A procedure of calibration of the micro-parameters is necessary, in which these are varied in order to produce the experimental values of the sample Young's modulus, Poisson's ratio and macroscopic strength. This fact introduces some difficulty in the application of BPM, as a set of preliminary runs is always needed, before a given practical problem is solved. This may be facilitated by automated procedures to replicate the lab test setups, available, for example, in PFC (Itasca, 2008). Neural network techniques have also been used to identify the micro-properties that reproduce the deformability and strength of rock samples in uniaxial compression tests (Tawadrous et al., 2009).
The analysis of the results of the simulation also requires tools unlike those in continuum modelling, for example, the definition of representative measures of strain and stress at various locations in the sample by means of measurement circles. Statistics of the evolution of performance indicators during the test, such as bond breakage in tension or shear, are also indispensable for a quantitative characterization of the response (e.g. Potyondy, 2012).

Modelling developments
The success of the BPM approach has promoted an active research effort in recent years, with many authors testing and enhancing the capabilities of particle models to reproduce more realistically the observed response of rock materials under the standard laboratory tests. The role of the numerical micro-parameters in modulating the various aspects of the physical response is now much better understood (e.g., Schöpfer et al., 2009).
In order to achieve a better match of the observed triaxial behaviour, Potyondy & Cundall (2004) had proposed that use of macro-particles composed of groups of circular particles joined together, therefore providing more interlocking in the assembly and enhanced frictional effects (Fig. 16). Subsequent work confirmed the success of this approach (Cho et al., 2007;Yoon et al., 2011;Potyondy, 2012). The ratio of uniaxial compressive strength to tensile strength in the numerical simulations, which was too low in the early models, also became easier to control. The macro-particles may behave as rigid bodies ("clumps"), or the bonds between the component particles may be allowed to break, therefore simulating the progression of intragranular cracks. These models essentially provided a way to discretize the grain structure in more detail, accomplishing more complex grain shapes and behaviour, naturally at the expense of increased computational effort.
Polygonal block models resemble more closely the grain structure of many rocks. However, from a computational point of view, they are much more demanding, mainly because the geometric contact calculations between the edges and vertices of polygons involve many more operations than those needed for circular shapes. Models with elastic blocks had already been used to study the effect of random joint patterns on stress distributions (e.g., Brady et al., 1986). Very interesting results are now being attained in fracture analysis. Damjanac et al. (2007) studied the micro-mechanical behaviour of lithophysal tuff specimens with both particles and blocks. Lan et al. (2010) represented the microstructure of brittle rock by means of a deformable polygonal grain-like assembly, in order to study effect of heterogeneity of grain deformability properties on the behaviour under uniaxial compression (Fig. 17). Kazerani & Zhao (2011) used both Voronoi and Delaunay block assemblies to study the fracture propagation in uniax-ial and Brazilian test specimens. Fracture of blocks is possible in some codes (e.g., Munjiza, 2004;Eberhardt et al., 2004), either through predefined potential crack paths, or by means of some form of internal remeshing. For this purpose, however, macro-particles formed by disks appear to be a more flexible and promising route, as they provide more freedom for the cracks to propagate with fewer numerical restraints.
The main disadvantage of polygonal models is the computational effort, as contact detection and particle interaction calculations are much slower than for circular particles. This fact limits the size of assemblies that can be handled with reasonable run times, particularly in 3D, making it difficult to approach scales larger than those of lab specimens. The proposal of formulations based on multiple point contacts between two circular particles provided a way to achieve a response closer to that of polygonal models in a much more efficient way (Azevedo & Lemos, 2005). In this approach the contact between two disks is assumed to take place along a line segment, and it is discretized into several points, which may fail independently, allowing progressive cracking to develop. A more gradual response is obtained, in comparison with the instant breakage of a single rotational bond, while the use of contact models with softening also provides a less brittle postpeak response. Similar concepts underlie the "flat-joint model" proposed by Potyondy (2012). A 3D version of the multiple contact model was presented by Azevedo & Lemos (2013), using concentric rings of contact points (Fig. 18). In practice, five contact points appear sufficient to give good results. A key aspect addressed in this paper is the importance of the assembly compactness, which may be measured by the coordination number, defined as the average number of interactions for each particle (in this respect a multiple point contact counts as one interaction). For polyhedral block assemblies this number is higher than for spheres. Therefore, the authors proposed a new criterion to detect contact between spheres. A Voronoi-Laguerre diagram was superimposed on the sphere assembly, and two particles were allowed to interact if the corresponding polyhedra were adjacent, even if the spheres were not exactly touching. This increased the number of interactions, providing a response closer to the one obtained with a much more expensive polyhedral block model, namely capturing the experimental failure envelope of Lac du Bonnet granite (Fig. 19).
The aim to create larger model sizes has led research into more efficient formulations, such as Cundall's (2011) "lattice model", in which the finite-sized particles are replaced by point masses, and the contacts between particles are replaced by springs that may break. Assuming small displacements, this model achieves high computational efficiency because the interaction geometry (location and apparent stiffness of springs) can be pre-computed, eliminating contact detection as an overhead.

Rock joint representation
In addition to rock fracture, other aspects of fundamental rock behaviour are being approached by means of detailed models at the meso-level scale, mostly using particle representations. The mechanical response of rock joints in laboratory tests is a suitable application for these models, capable of reproducing the roughness of joint walls. Joint shearing was simulated with PFC by Cundall (2000) in 2D, while a 3D model was used by Park & Song (2009). The complex patterns of compressive and tensile contact forces creating by the uneven contact conditions are evident in these analyses. The development of cracks in the vicinity of the joint walls during shearing of the samples may also be followed in these simulations (Fig. 20).
The effect of discontinuities on stress waves was investigated by Resende et al. (2011) also with a particle model. A simplified procedure was developed to create a rough joint, with uneven contact conditions, and a contact constitutive model with variable stiffness was employed, leading to a global stress-dependent behaviour. The propagation of a plane wave across the discontinuity for various in situ stress levels was then simulated. The differences in the contact response for the different normal stresses induced complex velocity fields in the vicinity of the discontinuity, creating marked differences in the reflected and transmitted waves obtained (Fig. 21).
The detailed modelling of rock joints has motivated much interest, namely to investigate the hydromechanical behaviour, for which various techniques were used to replicate the joint roughness and aperture distributions, generally assuming joint walls to be rigid or elastic. DEM particle models have the potential to improve the analysis of joint wall behaviour, addressing the effects of roughness wear and cracking. This is a promising research field, even if the current particle models still have considerable difficulty to attain a realistic resolution of roughness profiles, as well as the processes that lead to the progressive damage to joint surfaces observed during normal or shear loading tests.

Rock mass modelling
The main practical challenge posed to particle models today lies in their extension to problems of larger size, necessary to jump from lab test analysis to the scale of engineering structures. The question of its feasibility seems to depend mostly on forecasts of computing power increases. A broader question is whether this objective is truly necessary, or, more precisely, in which circumstances will these finer representations be superior to the existing simplified engineering models. Cundall (2001) Figure 18 -Multiple point contact for particle models (Azevedo & Lemos, 2013).

Figure 19
-Comparison of experimental and numerical envelopes in triaxial compression tests of Lac du Bonnet granite (after laws to simulate their complex behaviour, which tend to become excessively complicated, and involve many parameters with obscure physical meaning. Besides, cracking and failure in rock processes are more naturally approached by discontinuum methods. In contrast, Cundall argued, assemblies of bonded discrete particles are capable of capturing the complicated behaviour of actual material with simple assumptions and few parameters at the micro level, with complex overall behaviour arising as an emergent property of the assembly. There is no question that particle models are an indispensable research tool. However, in many practical engineering applications, simpler models appear sufficient for the purposes of safety assessment or prediction of global behaviour, as illustrated, for example, by the deformable block models presented above for dam foundation problems. After all, given the uncertainties involved in site characterisation, even these block models may challenge our capacity to provide reliable input parameters, so more refined representations may not be warranted.
One of the critical advantages of representing blocks as macro-particles composed of circular particles is the ease of considering block fracture, given the progresses discussed in the previous sections. The role of non-persistent discontinuities and the contribution of rock bridges to the strength of a rupture surface is an important aspect. Assumption of large extents for all joints is often overly conservative and contrary to geological information. The alternative to particle models is to adopt plasticity constitutive models for the block elements, which may be reasonable for global modes of response.
The Synthetic Rock Mass (SRM) concept proposes to extend the range of particle models to engineering scale, by considering the presence of the discontinuities (Mas Ivars et al., 2011). A discrete fracture network (DFN) is overlaid on a particle assembly, thus partitioning it into a system of grains or blocks formed by bonded circular particles, as shown in Fig. 22. Different properties are assigned to the bonds of the contacts between particles belonging to the same block, representing the intact rock material, and to the contacts between adjacent blocks, representing the joint behaviour. Joint may be persistent or terminate within blocks, thus allowing progressive crack development processes.  A major difficulty arises in the representation of sliding on the rough rock joints created in SRM models. The irregular shapes of the macro-particles, dependent on the size of the elementary disk or sphere components, imply excessive interlocking and dilation. From an engineering perspective, it would be desirable to assign a given friction angle, and controllable dilatant behaviour, to these surfaces. The solution to this problem lies in Cundall's Smooth Joint Model (SJM) concept, which is applied to the contacts between particles belonging to neighbouring blocks. Even if the interface is not an exact straight line, the SJM logic forces all these contacts to adopt a common normal, leading to a smooth sliding behaviour governed by the prescribed friction angle. Otherwise, sliding would imply riding over the particles, causing unrealistic dilation.
Successful applications of the SRM concept are starting to appear in the literature. Cundall (2007) presented a 2D model of a rock slope with a complex system of joint sets, which resulted in an assembly of about 330,000 particles. A detail of the model is shown in Fig. 23, displaying the contours of displacements triggered by the progressive excavation of the slope, which indicate the dominance of the toppling mechanism ensuing from the steeper dipping joints. A 3D SRM model based on the lattice formulation was applied by Cundall & Damjanac (2009) (Cundall, 2007).
employed it to study scale effects in jointed rock masses. The anisotropic response and the trends in tensile and compressive strength variation were investigated by performing a series of numerical tests on samples of various sizes. Pierce & Fairhurst (2011) discuss the use of SRM models in the investigation of caving processing in underground mines. The complex problem of hydraulic fracturing is also being approached with these fracture models (Damjanac et al., 2010). The expansion of particle models to larger scales, mostly limited by the available computing resources, is inscribed in a natural tendency towards more complex simulations. Its effective use, however, depends on our ability to fulfil some essential requirements. First, it is necessary to generate realistic macro-particle assemblies, in conformity to the specific joint patterns at the site. Robust analysis methodologies have to be established, namely involving expedite procedures for calibration of micro-parameters, and multiple simulation schemes to account for the variability of the physical data. Finally, the interpretation of the numerical results of complex models needs to be based on sound and objective criteria, appropriate to specific aims, for example, the safety assessment of foundations or excavations.

Closing Remarks
The widening range of problems facing rock engineering continues to stimulate the development of new analysis tools. Stricter safety and environmental requirements, as well as technological and economic factors, challenge our capability to predict the integrated behaviour of rock masses subject to the demands of new structures or excavations. Numerical models based on the discontinuum paradigm are one of the essential instruments available today. The development and application of discrete element models has followed two distinct but complementary paths. On one hand, we have the more sophisticated models employed in research, aimed at the understanding of the fundamentals of rock and rock mass behaviour, in conjunction with laboratory and well-controlled in situ experiments. On the other, the simpler but more robust techniques that suit the needs of engineering design practice. As the more complex models are progressively validated against experimental evidence, while becoming more reliable and userfriendly, they are increasingly called to address practical goals.
The effective use of advanced numerical tools relies, first of all, on our ability to supply geologic and geomechanical information to characterize adequately the site conditions. Statistical description of jointing and expedite methods for generation of numerical fracture networks are a key issue. The adequate representation of network connectivity, joint persistence, or rock bridges remains an important goal, in order to avoid overly conservative assumptions. Extending the numerical models to encompass multi-physics phenomena is perhaps the greatest challenge, given the range of problems involving the time-dependent interactions between hydromechanical, thermal, transport and chemical processes in rock masses. Discontinnum models are also a key tool to address the problems of hydraulic fracturing, dynamic rupture and fragmentation processes. Computational developments, namely in parallel processing technologies, appear vital to make large 3D simulations more manageable.
At present, engineering is faced with ever greater demands for a more rigorous evaluation and management of risk, viewed in its full human and economic aspects. The achievement of these tasks, however, rests on our ability to understand the behaviour of rock masses and engineering materials, and to predict their performance. Numerical modelling tools play a key role in this endeavour, always keeping in mind that the development of new techniques has to be complemented by attention to proper modelling methodologies and effective training in their application.