## Abstract

Naturally, deformed rocks commonly contain crack arrays (joints) forming patterns with systematic relationships to the large-scale deformation. Kinematically, joints can be mode-1, -2 or -3 or combinations of these, but there is no overarching theory for the development of the patterns. We develop a model motivated by dislocation pattern formation in metals. The problem is formulated in one dimension in terms of coupled reaction–diffusion equations, based on computer simulations of crack development in deformed granular media with cohesion. The cracks are treated as interacting defects, and the densities of defects diffuse through the rock mass. Of particular importance is the formation of cracks at high stresses associated with force-chain buckling and variants of this configuration; these cracks play the role of ‘inhibitors’ in reaction–diffusion relationships. Cracks forming at lower stresses act as relatively mobile defects. Patterns of localized deformation result from (i) competition between the growth of the density of ‘mobile’ defects and the inhibition of these defects by crack configurations forming at high stress and (ii) the diffusion of damage arising from these two populations each characterized by a different diffusion coefficient. The extension of this work to two and three dimensions is discussed.

## 1. Introduction

One of the outstanding problems in structural geology is the origin of joint (or fracture) patterns in deformed rocks. Commonly, joints occur in more than one set, each comprised multiple joints that share a common orientation, spacing and persistence with respect to other structural features, such as folds and faults (Norris 1967; figure 1). Understanding the mechanics that underlies the formation of such joint sets is of fundamental importance in areas such as hydrocarbon extraction, carbon sequestration and storage, geothermal energy extraction, underground radionuclide storage, contaminant transport in hydrological systems, the formation of ore deposits and the stability of open-pit and underground mines. Although the literature on joints is voluminous (Barton & Larsen 1985; La Pointe & Hudson 1985; Barton & Zoback 1990; Hobbs 1993), there is no underlying theoretical basis for understanding the patterns of joints that are commonly observed in naturally deformed rocks. In sedimentary and igneous rocks, joints are commonly mode-1 fractures that originate through shrinkage during compaction or thermal contraction during cooling. However, in deformed rocks, the joints can be mode-1, -2, -3 or a combination of modes, and several models for their formation have been proposed. These models range from the ‘single-event mechanism’ (Griffith 1920; Davidsen *et al*. 2007 and references therein), where one fracture propagates through the rock mass, governed by the laws of elastic fracture mechanics to models where damage accumulates locally comprised microfractures, and these damage concentrations finally coalesce to produce a throughgoing fracture (Hamiel *et al*. 2009). Models involving hybrids between these two extremes have also been proposed (Resches & Lockner 1994; Backers *et al*. 2005). All of these studies explore models for the initiation and final emergence of a single fracture. In this paper, we are concerned with a model that describes the growth of *patterns* in fracture arrays. There have been some studies of fracture pattern formation from a number of different viewpoints (Muhlhaus *et al*. 1996; Tang *et al*. 2008 and references therein), but, in general, an overarching theory is lacking.

We use the term *crack* to mean a fracture at a fine scale (normally the grain scale) that originates due to deformation. A crack is a form of *damage*, and damage accumulates in deformed rocks as a result of increasing crack densities. When the local damage density becomes large enough for individual cracks to coalesce to form a connected network, then the connected structure becomes a *joint*. The scientific issues include questions such as: ‘why is the network of cracks planar’ and ‘why do the cracks form patterns that ultimately lead to the growth of well-defined joint systems’? This paper is concerned solely with the second question.

In this paper, we explore the notion that joints in deformed rocks arise from the competition between the production of damage (by various mechanisms) and the diffusion of the density of this damage as the deformation proceeds. The approach is motivated by work on dislocation dynamics (Walgraef & Aifantis 1985*a*,*b*,*c*; Aifantis 1986; Zbib *et al*. 1996; Shizawa & Zbib 1999; Shizawa *et al*. 2001; Pontes *et al*. 2006). Such an approach suggests that the crack-pattern formation results from interactions during deformation between different crack populations, a process that tends to localize crack densities, and diffusion that tends to homogenize crack densities. Patterns develop when a critical stress is reached, corresponding to a balance between these two competing processes. The present work is based on damage evolution observed in numerical models comprised a large number of particles that are able to separate (mode-1 cracks), to shear (mode-2 cracks) and to interact to form more complicated arrays during deformation. The strain in these aggregates is accommodated solely by the formation of these cracks and diffusion of crack densities; the maximum density of cracks forms to accommodate the imposed strain rate. The study builds on a large number of previous studies employing distinct element approaches (see Cundall & Strack 1983; Tordesillas *et al*. 2004; Ord *et al*. 2007), in which patterning of damage evolution has been demonstrated. The traditional way of viewing these deformation patterns (figure 2*a*), motivated by the work on localization of deformation in rate-independent continua (Rudnicki & Rice 1975), is as an array of intersecting shear zones. The emphasis in this paper is on viewing the patterns as arrays of ‘crack cells’ by analogy with the work on ‘dislocation cells’ in metals (Holt 1970; figure 2*b*). The crack cell pattern is reminiscent of subgrains in deformed metals with rhomboid-shaped islands of less-deformed material containing only a few geometrically unnecessary dislocations separated by walls of tangled to highly-ordered dislocations. The dislocations comprising these walls are geometrically necessary to accommodate the deformation of the specimen. As with a dislocation microstructure, the patterns of normal and shear cracks have to diffuse in space and time to accommodate the deformation and deformation rate. They are likely also to diffuse at different rates. Understanding how crack densities are generated and diffuse in space and time is the essence of the problem.

In many papers on granular media (Kuhn 1999, 2005; Darve *et al*. 2004; Gardiner & Tordesillas 2004; Tordesillas *et al*. 2004), the distinct elements are circular (two-dimensional cylindrical) or spherical in shape, and in many there is no bonding between particles so that the material is cohesionless and behaves as an assembly of sand particles with only frictional interactions between the particles. We continue with circular (two-dimensional cylindrical) particles in this paper, but implement cohesion (both in shear and for normal displacements) as well as friction and elasticity between the particles (as in the early work by Cundall & Strack 1983). Independent of the presence or absence of cohesion between the particles, the same overall features arise during the deformation of these models:

Nonlinear elastic softening occurs prior to the development of plastic strain. This is associated with the development of microbands approximately parallel to the direction of shearing in inhomogeneous simple shearing experiments (Kuhn 2005) or at 40–60

^{°}to the maximum compression direction in a pure shearing deformation (Kuhn 1999). These represent local areas of dilation and compaction that develop from the onset of deformation, i.e. well within what one might commonly think of as the ‘elastic’ range of deformation. These structures are not persistent and tend to grow, interact and disappear. They appear to be unrelated to the strong force chains at any one time. Ord*et al*. (2007) show that the hysteresis in loading/unloading experiments for materials with cohesion is small in the nonlinear elastic regime, indicating that the formation of cracks is largely reversible, except that, although the cracks close upon unloading, the bonds between grains remain broken.Persistent shear bands tend to form relatively late in the deformation history (Kuhn 1999; Oda & Yoshida 1999). These are the classical shear bands that would be interpreted as the result of bifurcation in rate-independent continua (Rudnicki & Rice 1975; Vermeer & de Borst 1984). Oda & Konishi (1974

*a*,*b*), Cundall & Strack (1983), Oda & Kazama (1998), Walsh*et al*. (2007) and Hunt*et al*. (2010) have interpreted these as forming as a result of buckling or kinking of the strong force chains.There is commonly a large drop in stress (perhaps three-quarters of the total deviatoric stress) that accompanies the first development of localized deformation in the aggregate. This is marked by a localization of sites of rapid decrease in the potential energy of grains (Gardiner & Tordesillas 2004) and the first appearance of non-positive values of the second-order work (Darve

*et al*. 2004; Ord*et al*. 2007). There seems to be a consensus that localization coincides with the buckling of force chains in the aggregate (Cundall & Strack 1983; Ord*et al*. 2007; Walsh*et al*. 2007). The formation and evolution of strong force chains between grains (Cundall & Strack 1983) divides the system at a particular instant into ‘strong force-chain columns’ and ‘weak regions’ between these columns. Dissipation occurs when the strong force chains collapse, often arising from frictional sliding due to the motion of a single grain in a particular chain to occupy a neighbouring void (Tordesillas & Muthuswamy 2009). Such a collapse only occurs at high local stresses. Both translation and rotation of particles are involved in this process. However, this motion of a single grain influences the subsequent motion of grains up to 10 particle diameters removed from that grain.The zones of localized deformation are separated by an array of rhomboid-shaped regions (figure 2), where the strain is relatively small and the deformation is non-affine in the sense that plots of the departure from the homogeneous velocity field imposed on the aggregate show vortices. This non-affine deformation is a characteristic of many non-crystalline and disordered (Utter & Behringer 2008) materials, including glasses (Leonforte

*et al*. 2006; Del Gado*et al*. 2008). In the experiments reported here, these non-affine deformations contribute little to the total plastic strain, and seem to represent the coordinated shuffling of particles to accommodate heterogeneities in the deformation created by the localization of deformation into neighbouring zones. The cracks that form in the regions of non-affine strain tend to be mode-2 and hence involve frictional sliding.In materials where cohesion exists, the zones of localized deformation coincide with concentrations of cracks that are dominantly mode-1. Although breakage of bonds occurs between grains with subsequent frictional sliding, Wolf

*et al*. (1998) have shown that dissipation arising from frictional sliding occurs dominantly in the weak regions between strong force chains.

Over the past two decades, pattern formation has been investigated in a large number of physical, chemical, biological and geological systems (Nicolis & Baras 1984; Murray 1989; Ortoleva 1994; Walgraef 1997). One of the most successful ways of analysing and understanding such systems arises from the behaviour of various forms of coupled reaction–diffusion equations (Cross & Hohenberg 1993). In particular, Walgraef & Aifantis (1985*a*,*b*,*c*) and Aifantis (1986) have applied such concepts to the development of localized deformation patterns in metals. The results of these latter analyses indicate that such localization should develop, both as temporal oscillations in deformation intensity (a Hopf instability) and spatial fluctuations in deformation intensity (a Turing instability). Both of these kinds of instabilities are observed in deforming metals (Aifantis 1987).

The approach adopted here relies heavily on the above arguments developed by Walgraef & Aifantis (1985*a*,*b*,*c*, 1988) and by many subsequent workers, for the development of spatial and temporal patterning in dislocation densities in deformed metals. As such, we develop reaction–diffusion equations for the evolution of cracks in deforming particle assemblies with cohesion, and show that these lead to equations identical in nature to those of Walgraef & Aifantis (1985*a*) and to patterns of localized deformation characterized by spatial patterning of crack densities. In this paper, we first present some results of deforming an aggregate of particles in which there are elastic and cohesional inter-relations (both in shear and compression) between the grains, as well as frictional behaviour once the grains begin to slide. Emphasis is placed on the stress–strain response, the evolution of damage and the spatial and temporal patterns of localized deformation that form. We then develop a model for damage evolution based on these observations and the work on dislocation patterning in metals. Finally, we analyse this model in terms of a simple one-dimensional, reaction–diffusion equation, point out the implications for extension of the results to two and three dimensions and discuss the existence of size effects in geomaterials.

## 2. Deformation in granular aggregates with cohesion

The computational or virtual experiments described below are based on the physical experiments of Edmond & Paterson (1972). The material in this case uses the results for Gosford sandstone (see also Ord 1991; Ord *et al*. 1991, 2007; Hobbs *et al*. 2004) as a base-case scenario, and we then explore the effect of variations of the constitutive parameters used in that base case. The experimental work of Edmond & Paterson (1972) is used to constrain the constitutive parameters corresponding to the base-case computer simulations (Ord *et al*. 2007). The computations employ the distinct element code PFC2D (ITASCA 2004).

### (a) The virtual rock

The computer experiments reported here are two-dimensional and are performed in plane strain; the virtual material is 0.2 m square on edge. It is comprised particles with a range of radii from 2 to 3 mm, uniformly distributed and packed so as to approach a porosity of 0.13 (see figure 3) and a density of 2226 kg m^{−3}. A process analogous to physical compaction is followed (ITASCA 2004). Particles of artificially small radii are created within the specimen spaced at random coordinates. The particles are then expanded to radii within the desired range until the desired porosity is obtained, which, in this case, required 1772 particles.

For the base case (model 1), the initial normal and shear stiffnesses of the particle contacts, using a linear contact model, are both 1.16×10^{10} N m^{−1}, as is the wall normal stiffness, for wall–particle interactions, as described in Ord *et al*. (2007). Parallel bonds, providing the attribute of cohesion to the material, are created for all contacts. Both the normal and shear stiffnesses of the parallel bonds are 1.16×10^{10} Pa m^{−1}; and the parallel bond normal and shear strengths are 3.5 MPa. The frictional coefficient assigned to the particle surfaces varies between 0.5 and 1.0 (see Skinner (1969) for a discussion of the role of inter-particle friction in a spherical particle system). The material is pre-stressed to attain 10 MPa in both the vertical and horizontal directions. We have explored a range of parameters around the base-case scenario, and these values are given diagrammatically in figure 4 and in table 1. As indicated above, it has been found that the parameter suite used as the basis for model 1 reproduces the results derived by Edmond & Paterson (1972) from physical experiments on Gosford sandstone. We explore four other models in which the elastic modulus remains the same as for model 1, but the ratio of the normal-to-shear cohesive strength is varied as shown in figure 4 and table 1.

### (b) Boundary conditions

The two horizontal boundaries are rigid and are both moved inwards at a velocity of magnitude 1.6×10^{−2} m s^{−1}, while a confining pressure of 10 MPa is maintained on the vertical walls. The maximum total bulk shortening attained is 0.2. We report only results for strains up to 0.04 here.

The horizontal boundaries are accelerated gradually from zero velocity to the specified velocity before the formal experiment takes place; some loading of the specimen takes place at this stage. The loading rate during the formal experiment is sufficiently slow that the specimen remains in quasi-static equilibrium throughout the test (Ord *et al*. 2007). This was checked by monitoring the work done by the platens and by comparing this to the strain energy stored in the particle assembly (ITASCA 2004). These two energy quantities are required to be equal during the elastic portion of the experiment for the conditions to be quasi-static. Otherwise, the loading conditions are dynamic. It was found that this process of ensuring quasi-static deformation is very important. Quite different behaviours of force chains and patterns of particle motion develop if dynamic loading is imposed rather than quasi-static loading (Ord *et al*. 2007).

### (c) Mechanical behaviour

The bulk mechanical behaviour of the virtual specimen is represented by the stress–strain curves shown in figures 5 and 6 for model 4 (see tables 1 and 2) and in figure 7 (right-hand side column). The beginning of the formal experiment is assumed to be that step when the velocity specified for the loading platens for the experiment is attained, and this is represented by zero axial strain. The peak or maximum deviatoric stress (table 2) attained by the virtual specimens ranges from about 21 MPa for model 1 to 18 MPa for model 5.

The non-affine velocity field and the continuous lines of contact forces (force chains) for model 4 are shown at axial strains of 0.008, 0.02 and 0.04 in figure 5. Similarly, figure 6 shows the progressive development of localized zones of shearing and fracturing throughout the specimen, as well as the evolution of the number of normal and shear cracks (figure 6*h*). These results may be compared with those for model 1 that are given by figure 3 of Ord *et al*. (2007). Figure 7 shows the different fracturing and mechanical behaviours associated with the five different models represented by figure 4.

The unloading behaviour of the material represented by model 1 is represented by fig. 4 in Ord *et al*. (2007). There is minor hysteresis of the unloading/loading curve at low strains. The elastic modulus stays very much the same. This is followed by a nonlinear response with little hysteresis until after the elastic modulus of the material begins to soften at a strain of about 0.008. Thereafter, hysteresis is marked, and the gradient of the unloading curve decreases with increasing strain. Again, the volume strain/axial strain behaviour (Ord *et al*. 2007, fig. 4*b*) shows convincingly the elastic response at low strains, and strong hysteresis in the volume strain at higher strains.

The microstructural responses of the models are illustrated in figure 7 for models 1–5 and for axial strains of 0.02 (second column) and 0.04 (third column). The distribution of normal (mode-1) and shear (mode-2) cracks is shown for each model. In model 1, where the constitutive parameters are identical for both normal and shear failures, normal cracks dominate and are mimicked by the distribution of shear cracks. Slight changes in these constitutive parameters (compare models 4 and 5) result in significant differences in microstructural behaviour. In model 4, where the shear cohesion parameter is reduced, the microstructural behaviour is similar to model 1, but with a more diffuse distribution of normal and a greater density of shear cracks. In model 5, where the normal cohesion parameter is reduced, normal cracks dominate with very few shear cracks forming and no localization. Large changes in cohesional properties with respect to model 1 result in a dominance of one mode of cracking or the other with no localization (at least up to a strain of 0.04). Thus, for model 3, where the normal cracks are weakest, normal cracks are the only cracks formed and there is no localization. For model 2, the shear cracks are weakest and dominate with few normal cracks by a strain of 0.04 and weak localization.

Figure 7 shows that localization is favoured if normal (mode-1) cracks dominate, and this is the case if the densities of both normal and shear (mode-2) cracks are within an order of magnitude of each other (models 1 and 4). If only normal cracks develop (model 3), then no localization occurs; this corresponds to the normal cohesive strength being much less than the shear cohesive strength. This is also true for model 5, where normal cracks dominate, but the normal and shear cohesive strengths are closer together. If shear cracks dominate with very few normal cracks (model 2), which corresponds to the shear cohesive strength being much less than the normal cohesive strength, then localization is absent or diffuse. Localization depends on the development of large densities of both normal and shear cracks (models 1 and 4). This situation supplies opportunities for interaction and the formation of more complex arrays, as illustrated in figure 8.

## 3. A model for the development of crack patterns

Motivated by the work of Walgraef & Aifantis (1985*a*,*b*,*c*), we propose that the formation of spatial patterns in the deformed aggregates described above can be analysed in terms of the evolution of two or more populations of cracks, as described below. As far as possible, the notation of Walgraef & Aifantis (1985*a*,*b*,*c*) is used here, although the terms have meanings that are different in most cases. The production of crack populations is treated as though these processes are chemical reactions. This means that cracks develop from the uncracked state and/or react with other cracks to form new populations of cracks; these crack densities can diffuse with diffusivities based on the configuration of the crack population. Single mode-1 cracks can diffuse with relative ease since such a process involves only the breakage of a neighbouring bond and can be accommodated in the surrounding network by elastic distortions. The diffusion of mode-2 cracks is intrinsically more difficult and involves not only the breakage of a neighbouring bond, but frictional slip on adjacent newly formed cracks. The diffusion of more complicated arrays, such as those generated by force-chain buckling, involves plastic distortion of neighbouring groups of grains and hence is more difficult again. We elaborate on these processes later in §3.

If we treat the generation and interaction of cracks as being formally equivalent to chemical reactions, then the thermodynamic treatment of the coupling between deformation, chemical reactions and thermal transport, as developed by workers such as Coussy (1995, 2004), is relevant. We outline this approach below and discuss a general non-isothermal deformation with coupled reactive, diffusive and thermal behaviour. This is done to set the scene for a more general consideration of what may occur in natural deformations in §4. We concentrate on isothermal deformations with respect to the granular aggregates considered here. Symbols used in the following discussion are summarized in table 3, together with their units.

In discussing the thermodynamics of chemically reacting systems where the reaction has not proceeded to completion, it is convenient to use the state variable, *ξ*, called the extent of the reaction with 0≤*ξ*≤1 and the conjugate variable, *A*, the affinity of, or driving force for, the chemical reaction (Kondepudi & Prigogine 1998). is the rate of the reaction, where the overdot indicates the material derivative with respect to time *t*. *A* is a linear function of the difference in the sum of the chemical potentials of the reactants and products in the reaction (Kondepudi & Prigogine 1998). The extent of a reaction can also be measured by the concentration of product development during the reaction or by the depletion in concentration of a reactant. In the context of crack patterning, the ‘reaction’ corresponds to crack generation or crack depletion in an interaction that produces a more, or less, complicated array. In addition to the extent of the reaction, one can use the density, *ρ*_{K}, of the particular crack population, *K*, as a variable that defines the thermodynamic state of the system.

Non-equilibrium thermodynamics attempts to describe the evolution of systems not at equilibrium in terms of two potentials, the Helmholtz free energy (or if pertinent to the problem, the Gibbs free energy) and the dissipation function (Houlsby & Puzrin 2006). For instance, if we consider the specific Helmholtz free energy, *Ψ*, as the relevant function, then we can choose the thermodynamic state variables as the elastic strain, , the absolute temperature, *T*, the extent, *ξ*_{K}, of a reaction that produces a given crack population, *K*, and the density, *ρ*_{K}, of the *K*th crack population (Kondepudi & Prigogine 1998). We write
3.1
where *u* is the specific internal energy and *s* is the specific entropy.

From equation (3.1) and following standard arguments in continuum thermodynamics (Callen 1960; Coussy 2004), we obtain
3.2
where *σ*_{ij} is the Cauchy stress and *μ*_{K} is the ‘chemical potential’ of the *K*th component. Here, the ‘chemical potential’ is interpreted as the mechanical energy required to produce a crack, or a crack array, in the stressed solid. *v*_{0} is the specific density of the material and is the initial density of *K* cracks.

To proceed, we need to couple temperature to these processes and derive an expression (the energy equation) that describes the ways in which the various dissipative processes contribute to temperature changes. A combination of the first and second laws of thermodynamics may then be written (Coussy 2004) as
3.3
where Φ is the dissipation function, and Φ^{mechanical}, Φ^{diffusive}, Φ^{reaction} and Φ^{thermal} are the contributions to the total dissipation from purely mechanical processes, diffusion, crack reactions and thermal diffusion, respectively. The mechanical dissipation includes the energy required to introduce a density, *ρ*_{K}, of the *K*th crack population into the deforming material through a crack reaction, whereas the diffusive term includes the flux, *J*_{K} of the *K*th crack species across gradients in the chemical potentials, *μ*_{K}, and the temperature.

These dissipation functions are given by (Coussy 2004)
3.4a
3.4b
3.4c
3.4d
where summation on *K* is not intended, but *K* takes on the values 1,…,ℵ. If we intend summation on *K*, then we make it explicit by use of the S symbol. In equation (3.4b), the mid-dot represents the scalar product of vectors; in equation (3.4d), *κ*^{thermal} is the thermal diffusivity and *c*_{p} is the specific heat at constant pressure. If equation (3.4c) was to express a chemical reaction, then another term involving the rate of heat generated by that reaction would need to be added.

Using the expressions *s*=−∂*Ψ*/∂*T* and , equation (3.4a) becomes
3.5
If we assume, ignoring any possible effect of thermal expansion/contraction on the generation and propagation of the cracks, that
3.6
then using equation (3.3), we can eventually arrive at the energy equation that expresses the change in temperature arising from all of the dissipative processes and hence does not include contributions from elastic deformations,
3.7
where *χ* is the Taylor–Quinney coefficient and represents the proportion of mechanical work arising from dissipative deformation that is available to increase the temperature or to drive diffusion and crack reactions. At high strains, where the energy arising from deformation is stored in crystal defects, *χ* is generally in the range 0.85≤*χ*≤1 (Taylor & Quinney 1934).

Equation (3.7) is quite general and has been included as it points the way to a unified theory of pattern formation in joints formed in natural deformations where thermally activated processes and true chemical reactions operate synchronously with cracking. We return to this in §4. For the purposes of discussion of the granular aggregates used here, however, the deformation is isothermal and the only contribution to plastic deformation is the breaking of bonds between particles, the subsequent sliding between particles and the diffusion of damage. Thus, equation (3.7) reduces to ℵ coupled reaction–diffusion equations (Cross & Hohenberg 1993) of the form 3.8

In the context of crack generation, and are to be equated with the density rate and the reaction rate that generates the *K*th group of cracks, respectively, so that (3.8) is of the form
3.9

We distinguish two families of cracks: those that can form relatively easily and those that need a higher stress or strain to form. The ease of formation, as indicated in figure 7, depends on the bond strengths (shear and normal). Thus, those with low bond strength form more easily than those with high bond strength. These cracks are not restricted to be simply mode-1 or mode-2 cracks, and may be combinations of both. For convenience and by analogy with the dislocation literature, we call those that are easy to form, mobile cracks, and those that require a higher stress to form, immobile cracks. Each of these families are generated by different mechanisms, but once formed can diffuse through the matrix with different diffusion coefficients. They can also interact to form more complicated arrays so that a mobile crack can be ‘pinned’ within an array of immobile cracks, and an immobile crack can be liberated from an array to become a mobile crack. Thus, for each family, there is a balance equation of the form (3.9) that relates the generation, annihilation, interaction and diffusion of crack densities for that family and the relation of these processes to the evolution of the crack densities for other families. These processes and interactions are summarized below.

*Diffusion*of cracks takes place by what would be called in the metals literature, a stress-assisted vacancy diffusion mechanism. Thus, a mode-1 crack can diffuse to a neighbouring site by motion of a single particle, as shown in figure 8*a*. This involves breaking a bond, but can be accommodated in the network by elastic distortion of the neighbouring particles or by the creation of other cracks between adjacent particles. At low stresses, this diffusive motion can be achieved by elastic distortions, but at high stresses, the process is stress-dependent since it involves the breakage of bonds between grains. This process gives rise to a term*D*_{K}(∂^{2}*ρ*_{K}/∂*x*^{2}) in the balance equation for the family*K*.*D*_{K}is the diffusion coefficient for that family and*x*is a spatial coordinate. In what follows,*K*takes the identity of either mobile cracks, ‘m’, or immobile, ‘im’, cracks.*Pinning*of cracks takes place when large displacements of individual particles occur. The microstructures illustrated in figure 7 indicate that the type of transition illustrated in figure 8*b*is common. This only occurs at high stresses, and consists of a single particle moving sideways to generate two mode-2 cracks and a mode-1 crack. This, in fact, is the transition proposed by Tordesillas & Muthuswamy (2009) for buckling of force chains. If the mode-2 cracks have mechanical properties such that they only form at high stress (which is the situation for buckling of force chains), while the mode-1 cracks form at lower stresses, then this transition is equivalent to adding a sink term to the balance equation for mobile cracks and a source term to the balance equation for immobile cracks.*γ*is the rate at which these interactions occur. Other more complicated interactions are indicated by the observed microstructures (figure 7) that could add terms of higher order. Schiller & Walgraef (1988) show that the introduction of such higher order terms makes little difference to the final outcomes.*Liberation*of cracks occurs when an array, such as that formed by force-chain buckling, is relaxed by the motion of adjacent particles, such that an immobile crack can be converted to a mobile crack and diffuse away from the relaxed array, as shown in figure 8*c*. This is equivalent to adding a source term +*b**ρ*_{m}to the balance equation for mobile dislocations and a sink term −*b**ρ*_{m}to the balance equation for immobile cracks.*b*is the rate at which such a transition occurs. Notice that*b*is a parameter that only switches on at high stress, so that it has small values (or zero) at low stresses and only switches on when a critical stress is attained, which initiates the generation of complicated arrays of cracks such as occurs in force-chain buckling. By analogy with dislocation dynamics (Schiller & Walgraef 1988), an upper value of*b*is set by the imposed strain rate that can only be achieved by liberation of mobile dislocations. Some of the crack reactions are such that it is very difficult for some of the cracks in the interaction array to be liberated by diffusion. These cracks presumably remain pinned throughout the subsequent deformation.Initial

*generation*of cracks occurs at low stresses by the*formation*of mode-1 and mode-2 cracks in response to non-affine deformations (figure 8*e*). This damage is equivalent to ‘statistically stored’ dislocations, as initially proposed by Ashby (1969) for dislocations in the non-homogeneous plasticity of materials. Such cracks are randomly oriented and are the response to shuffling of individual grains to accommodate local inhomogeneities in the imposed deformation. Early in the deformation, these accommodate mainly nonlinear elastic deformation so that much of the crack opening is reversible (see fig. 4 of Ord*et al*. 2007), although no healing is incorporated in these models, so they remain as (closed) cracks. At later stages in the deformation, mode-1 and -2 cracks, constituting geometrically necessary cracks, form to accommodate the imposed deformation (figure 8*d*,*f*). The generation of these cracks is associated with plastic yield of the material and is expressed by a function*g*(*ρ*_{im}). The mathematical form of*g*is discussed by Schiller & Walgraef (1988). We propose, in keeping with Schiller and Walgraef, that the simplest form of the term is , where*ρ*^{0}_{im}is the initial density of immobile cracks and*a*is the rate at which immobile cracks are generated. This proposes that immobile cracks are initially generated in proportion to the initial density of defects and*g*takes on a positive value, but as the deformation continues, this rate is offset by the annihilation of arrays, so that the rate of generation of immobile cracks becomes negative, thus imposing an upper limit to the crack density.

The above discussion indicates that the overall evolution of these crack arrays is described by two coupled reaction–diffusion equations that express the relationships and coupling between diffusion, generation, annihilation and pinning of cracks,
3.10a
and
3.10b
This is the standard set of reaction–diffusion equations introduced by Walgraef & Aifantis (1985*a*), and has been intensively studied over the past two decades. In particular, this coupled set of equations has the following properties:

the homogeneous, steady state of the system is defined by 3.11

instability, represented by temporal oscillations (a Hopf instability) of the crack densities, occurs when 3.12

an instability represented by spatial patterning (a Turing instability) occurs when the stress becomes high enough that 3.13

the wavevector for such spatial patterning is 3.14 where

*λ*_{Turing}is the wavelength of the patterning; andthe Turing instability is reached before the Hopf instability if 3.15

This is true if *D*_{im}≪*D*_{m} and so is always true in these models.

Thus, the analysis indicates that an initial homogeneous distribution of cracks will spontaneously develop into a spatial pattern whose wavelength depends on the diffusion coefficients of the two families of cracks, the rates of generation (*a*) and of interaction (*γ*) of cracks and the initial density of immobile cracks (*ρ*^{0}_{im}). If the mechanical properties of bonds is such that only one family of cracks is generated, then no patterning is expected, as is shown in figure 7 (models 3 and 5).

## 4. Implications of the model

The spacing of localized zones of deformation in the model described here is controlled by the parameters *a*, *γ*, *ρ*^{0}_{im}, *D*_{im} and *D*_{m}. Even for dislocation arrays where a dynamical theory of dislocations exists, it is difficult to arrive at rigorous estimates of these quantities. For crack formation where no such theory exists, it is even more difficult to be rigorous in deriving values for these parameters. Some estimates are made below.

We consider model 4 shown in figure 6, which shows that the density of cracks diffuses over approximately 0.02 m in 0.125 s, giving a diffusivity of 3.2×10^{−3} m^{2} s^{−1}. We assume that this gives the value of *D*_{m}, and we take *D*_{im} to be two orders of magnitude smaller to be consistent with the theory giving *D*_{im}=3.2×10^{−5} m^{2} s^{−1}. Following Schiller & Walgraef (1988), we take the fraction *a*/*D*_{im} to be the inverse square of the mean free path of immobile cracks, and we take this to be equal to the grain size (3×10^{−3} m). Thus, *a*=3.6 s^{−1}.

At the stage where the stress–strain curve first starts to be nonlinear (∼0.005 strain), the total crack density is about 2.5×10^{2} m^{−2}. Again, we assume this is the value of and assume m^{−2}. If we assume that *b* has a value less than *b*_{Hopf}, an upper limit for *γ* is given by combining equations (3.11*b*) and (3.12) to give *γ*=6×10^{−3} s^{−1}. These values give an upper limit for the value for *λ* of 0.19 m. This compares favourably with observations, the shear zones being spaced at about 0.1 m. Moreover, this gives a value of *b* of 3.6 s^{−1}. Equation (3.13) also gives *b*=3.6 s^{−1}, so these estimates are consistent, and indicate that a value for *λ* of 0.19 m is perhaps reasonable.

This discussion has been concerned with one-dimensional arguments, but the extension to two and three dimensions is considered by Walgraef & Aifantis (1985*a*) and Aifantis (1986). The basic principles discussed here hold for these cases, and, in particular in two dimensions, the spacing and degree of development of localized zones of damage are expected to depend on the degree of anisotropy in the diffusion coefficients *D*_{m} and *D*_{im} (Aifantis 1986). Thus, whether one set of joints is developed or not in deformed rocks could be a reflection of this anisotropy.

In natural deformations, a number of processes other than isothermal crack formation are evident. Thermal effects become important due to the production of heat during the generation and sliding of cracks (Wu *et al*. 2002). The coefficient of friction and the cohesion may be temperature-dependent through an Arrhenius relation (Veveakis *et al*. in press). In addition, chemical reactions occur during the formation of cracks or during subsequent sliding. This includes the absorption of OH onto newly deformed surfaces (Zhu *et al*. 2005), and the precipitation of quartz and hydrous minerals such as chlorite from fluids circulating within the cracks. Some of these chemical reactions are endothermic (precipitation of quartz) and others are exothermic (OH absorption and precipitation of hydrous minerals), so that the complete set of equations labelled (3.4a–*d*) become relevant, and it is understood that the ‘reactions’ considered in (3.4a–*d*) involve not only the damage reactions discussed above, but strict chemical reactions also.

The equations that result from such considerations are reaction–diffusion equations, and some progress in considering these has been made by Veveakis *et al*. (in press). It turns out that equations of this type have been widely considered in the literature on combustion physics (Law 2006). The important point is that these approaches show promise, not only in understanding aspects of instability during frictional sliding (Hopf instabilities), but also in developing an understanding of the patterns developed through fracturing in naturally deformed rocks (Turing instabilities).

One aspect of natural fracture patterns is the proposal that these patterns have a fractal geometry. The theory presented here gives a scaling law, which results directly from (3.14), indicating that the wavelength of the fracture patterns is inversely proportional to the square of a length scale. For a true fractal relationship to exist, this length scale must also be a function of other length scales. For instance, Muhlhaus *et al*. (1996) show that a uniform array of parallel wing cracks in compression becomes unstable at a spacing that is a nonlinear function of the length of the mode-1 cracks in the wing crack and of the coefficient of friction. A bifurcation occurs at a critical value of the ratio of wing-crack spacing to wing-crack length so that a periodicity develops in the spacing of cracks. This process is driven by the elastic interactions between cracks, and it adds another length scale into the evolution of the crack array. These kinds of processes that depend on crack interaction in nonlinear ways presumably combine with expressions such as equation (3.14) to produce a fractal spatial pattern of cracks.

Lastly, Glazov & Laird (1985) have indicated that the theory, as developed by Walgraef & Aifantis (1985*a*), is strictly applicable to an infinite body. If the body is of finite size characterized by a length scale *L*, then equation (3.14) needs to be modified by a factor that is proportional to *L*. This effect may go some way to understanding the size effect seen in geotechnical applications (see da Cunha 1990 for references).

## 5. Concluding remarks

We have presented a model for pattern development of cracks in deforming granular media with cohesion, based on work in metals by Walgraef & Aifantis (1985*a*) and on the microstructures observed in deforming granular media. The crack patterns develop due to competition between the formation, interaction and diffusion of ‘mobile’ and relatively ‘immobile’ cracks or crack arrays. In particular, the immobile crack arrays are those associated with force-chain buckling, which is commonly attributed to localization phenomena in granular media. The patterns that develop are viewed as Turing instabilities, which develop once a critical stress is reached that enables complicated crack interactions to develop. The wavelength of the patterns predicted by the theory agrees with observations. Moreover, patterns do not develop unless the constitutive parameters for the material allow more than one population of cracks to form.

The model has a basis in non-equilibrium thermodynamics that enables an extension to natural fracture patterns which also involve thermal dissipation and endothermic/exothermic chemical reactions during their development.

## Footnotes

One contribution of 17 to a Theme Issue ‘Patterns in our planet: applications of multi-scale non-equilibrium thermodynamics to Earth-system science’.

- © 2010 The Royal Society