EGU General Assembly 2020
© Author(s) 2020. This work is distributed under
the Creative Commons Attribution 4.0 License.

A small-scale numerical study of fault slip mechanisms using DEM

Nathalie Casas1,2, Guilhem Mollon1,3, and Ali Daouadji2
Nathalie Casas et al.
  • 1LaMCoS, INSA Lyon/University of Lyon/CNRS UMR5259, Bâtiment Sophie Germain, 27bis Avenue Jean Capelle, 69621 Villeurbanne Cedex, France
  • 2GEOMAS,INSA de Lyon/University of Lyon, Bâtiment J.C.A. Coulomb, 34 avenue des Arts, 69621 Villeurbanne Cedex, France
  • 3Laboratoire de Géologie, École Normale Supérieure/CNRS UMR 8538, PSL Research University, 24 rue Lhomond, F-75005 Paris, France

How do earthquakes start? What are the parameters influencing fault evolutions? What are the local parameters controlling the seismic or aseismic character of slip?

To predict the dynamic behaviour of faults, it is important to understand slip mechanisms and their source. Lab or in-situ experiments can be very helpful, but tribological experience has shown that it is complicated to install local sensors inside a mechanical contact, and that they could disturb the behaviour of the sheared medium. Even with technical improvements on lab tools, some interesting data regarding gouge kinematics and rheology remains very difficult or impossible to obtain. Numerical modelling seems to be another way of understanding physics of earthquakes.

Fault zone usually present a granular gouge, coming from the wear material of previous slips. That is why, in this study, we present a numerical model to observe the evolution and behaviours of fault gouges. We chose to focus on physics of contacts inside a granular gouge at a millimetre-scale, studying contact interactions and friction coefficient between the different bodies. In order to get access to this kind of information, we implement a 2D granular fault gouge with Discrete Element Modelling in the software MELODY (Mollon, 2016). The gouge model involves two rough surfaces representing the rock walls separated by the granular gouge.

One of the interests of this code is its ability to represent realistic non-circular grain shapes with a Fourier-Voronoï method (Mollon et al., 2012). As most of the simulations reported in the literature use circular (2D)/spherical (3D) grains, we wanted to analyse numerically the contribution of angular grains. We confirm that they lead to higher friction coefficients and different global behaviours (Mair et al., 2002), (Guo et al., 2004).

In a first model, we investigate dry contacts to spotlight the influence of inter-particular cohesion and small particles on slip behaviour and static friction. A second model is carried out to observe aseismic and seismic slips occurring within the gouge. As stability depends on the interplay between the peak of static friction and the stiffness of the surrounding medium, the model includes the stiffness of the loading apparatus on the rock walls.

The work presented here focuses on millimetre-scale phenomena, but the employed model cannot be extended to the scale of the entire fault network, for computational cost reasons. It is expected, however, that it will lead to a better understanding of local behaviours that may be injected as simplified interface laws in larger-scale simulations.

How to cite: Casas, N., Mollon, G., and Daouadji, A.: A small-scale numerical study of fault slip mechanisms using DEM, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-8976,, 2020

Comments on the presentation

AC: Author Comment | CC: Community Comment | Report abuse

Presentation version 1 – uploaded on 29 Apr 2020
  • CC1: Comment on EGU2020-8976, Paul Selvadurai, 03 May 2020

    Thank you for this interesting prestenation.

    I was wondering about the contact laws (slides 15 and 16). You state: "Unbroken bond: constant value of cohesion - Broken bond: only inter-particular friction (µmicro=0.5) At the beginning of the experiment, all the particles in contact receive a percentage of cohesion. Once broken, a cohesive contact cannot be cohesive anymore"

    Q1: Can you relate this to a physical process observerd in the field of contact mechanics or tribology? Could this be considered a type of wear analogy?

    • AC1: Reply to CC1, Guilhem Mollon, 04 May 2020

      Dear Paul, thank you for your question. Here the irreversibility of the contact law (i.e. the fact that it goes from an "intact" to a "broken" status with a definitive loss of cohesion) is meant to represent a local damage of the cemented bonds between gouge grains. It leads to a weakening of the fault with slip, and will enable us to trigger dynamic instabilities in future versions of the model.

      It indeed bears some simillarities with tribological models, since we are part of a tribology lab. We commonly use such models in industrial applications to represente the emission of "third body" (i.e. damaged or worn material) in a contact interface, and the subsequent changes in the mechanical behavior of the contact.

  • CC2: Comment on EGU2020-8976, Jean Paul Ampuero, 04 May 2020

    Would you mind expanding on the "cohesion quantification" introduced in slide 20 (or point me to a reference explaining it in more detail)? I am not sure I understand the definition. Thank you.
    JP Ampuero

    • AC2: Reply to CC2, Nathalie CASAS, 04 May 2020

      Hi again, and thanks for your question:

      We have taken a "reference" or "limit value" of the rock we consider with the higher energy release. (here the Chilhowee quartzite) From the energy and the contact laws of our model, we extracted a value of cohesion. And we say that this cohesion is the maximum value we can found in real rocks, as if all the grains would be surrounded with cohesion. So we consider this value of cohesion as the maximum one. The cohesion in every simulation is then expressed in a percentage of this maximale cohesion. (We have done that because, the cohesion we enter in the simulation is a numerical one, and we wanted to relate it to real values )


    • AC3: Reply to CC2, Guilhem Mollon, 04 May 2020

      Dear Jean-Paul, thank you for your question. Each cohesive bond that we set up between any pair of contacting grains in the initial state of our model requires a certain amount of mechanical energy for breaking (related both to the bond tensile/tangential stiffness and to its tensile/tangential strength).

      To quantify the total amount of cohesion energy in the initial state of the fault (i.e. the energy that would be needed in order to break all the bonds), we normalize it with respect to what we consder as a physical upper limit for such a total energy. This upper limit corresponds to a surface energy of 62J/m² applied on the whole external surface of all the grains present in the simulation. This would correspond to 100% of cohesion, and we explore the influence of this percentage. Do not hesitate to ask if I was not clear enough.

      • CC5: Reply to AC3, Jean Paul Ampuero, 04 May 2020

        Than you both. I think I understand the definition now, in particular that the reference energy is integrated across all the grain boundaries. This motivates a follow-up question: In reality the grains are not in cohesive contact over their whole perimeter. So, in your simulations, what is the typical ratio of inter-grain contact surface to total grain surface? 

        • AC7: Reply to CC5, Guilhem Mollon, 04 May 2020

          This is a good point. The very notion of contact surface (or length, in 2D) is dubious in DEM since we are more talking of contact "points". Contacts are handled on a node-to-segment basis. To simplify, we consider that the length of a given contact between two grains is equal to the length of the contacted segment.

          A typical grain has 20-40 such segments (Nathalie, correct me if I'm wrong), and it also has 5-8 contacts with neighbouring grains (depending on the compation done prior to initialization of cohesive bonds). Hence, we could say that, in the initial state, 10-30% of a grain contour is in a bonded contact with its neighbours. This is meant to represent cemented bridges with a certain spatial extension.

          At the end of the day, simulations seem to indicate that what really matters is the energy needed to break each bond, not its "length". Which is why we are using this metric to quantify the initial state of cohesion of the fault.

          • AC8: Reply to AC7, Nathalie CASAS, 04 May 2020

            Yes closer to 40 segments !

  • CC3: Comment on EGU2020-8976, Carsten Uphoff, 04 May 2020


    how do you handle contact detection efficiently given that your grains are non-convex?


    • AC4: Reply to CC3, Guilhem Mollon, 04 May 2020

      Dear Carsten,

      the contact detection is performed on a node-to-segment basis. Each "slave" node of a given grain detects the proximity (and possible penetration) of each of the "master" segments of a neighbouring grain. This is common in FEM-based contact mechanics, but rather uncommon in DEM. A two-pass approach is used, such that each grain is both "slave" and "master". This way, grains with arbitrary shapes can be used.

    • AC5: Reply to CC3, Nathalie CASAS, 04 May 2020

      Hi, thank you for your question. We use a node-to-segment formulation based on three stages of detection : a broad proximity detection (to define pairs of close grains), a close proximity detection (to find contacts between a node from one particle with the segment of an other one), these two first steps are not performed at every time step, we have to choose the time period in order to make sure that we can not miss a contact and to minimize computational cost !The last step plays at every time step ( to perfome the actual contact detection)--> for more details here

      Guilhem, do you want to add something ?

      • CC4: Reply to AC5, Carsten Uphoff, 04 May 2020

        Thank you for your answers and the link to the paper. But basically if you have a polytope with N segments then the complexity of close proximity detection (Fig. 4) is O(N^2) as you need to test every vertex with every segment, right?

        I was wondering whether convex decomposition of your grains is an option such that you can use fast algorithms such as the GJK algorithm (

        • AC6: Reply to CC4, Guilhem Mollon, 04 May 2020

          Indeed, the complexity of this step is N². It is not a problem for us usually because N remains small (never more than a few dozens points on the contour of a grain), and because we only do it periodically in the simulation (like every 10-100 computation steps) to refresh the proximity lists of each object. Also, it parallelizes very well since it is a grain-wise.

          At the moment we have never needed to use lower-level convex objects. Actually, our contact algorithm was designed for highly-deformable grains (where contact handling has a negligible cost compared to the rest), and more effort was put on robustness than on efficiency.

          When going to 3D, in future versions of the code, it could be useful to consider your suggestion more closely, as detection costs might explode. We are considering HPC but are not yet skilled for this computation scale. Thanks for suggesting.