10626


FACTA UNIVERSITATIS  

Series: Mechanical Engineering  

https://doi.org/10.22190/FUME220327020F 

© 2020 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND 

Original scientific paper 

NUMERICAL SIMULATION OF DYNAMICS OF BLOCK MEDIA 

BY MOVABLE LATTICE AND MOVABLE AUTOMATA 

METHODS 

Alexander E. Filippov1,2, Valentin L. Popov2 

1Donetsk Institute for Physics and Engineering, Donetsk, Ukraine 
2Technische Universität Berlin, Berlin, Germany  

Abstract. Two versions of modified Burridge-Knopoff model including state dependent 

friction, elastic force and thermal conductivity are derived. The friction model 

describes a velocity weakening of friction and elasticity between moving blocks and an 

increase of both static friction and rigidity during stick periods as well their weakening 

during motion. It provides a simplified but qualitatively correct behavior including the 

transition from smooth sliding to stick-slip behavior, which is often observed in various 

tribological and tectonic systems. Attractor properties of the model dynamics is studied 

also. The alternative versions of the model are proposed which apply a simulation of 

the motion of interacting elastically connected mesh elements and motion of relatively 

large solid blocks, utilizing technique of the movable cellular automata. First version 

of the model was already basically studied before. Its advanced version here involves 

all components of the real system: state-depending friction and changeable rigidity, as 

well as heat production and thermal conductivity. Model based on the movable 

automata also involves the components included into traditional lattice model. It has 

its own ad-vantages and disadvantages which are also discussed in the paper. 

Key words: Burridge-Knopoff model, Friction, Tectonics, Voronoi diagram, Delaunay 

triangulation 

1. INTRODUCTION  

One of the most important applications in the study of the behavioral patterns of 

mechanical systems of block structures are the geological faults and block environments. 

Of particular interest, of course, is the study of the possible variations of their 

deformation, potential prediction or even prevention of the strong seismic shocks. The 

                                                           
Received: March 27, 2022 / Accepted April 23, 2022 

Corresponding author: Alexander E. Filippov 

Donetsk Institute for Physics and Engineering, R. Luxemburg Str. 72, 83114 Donetsk, Ukraine 

E-mail: filippov_ae@yahoo.com 



2 A.E. FILIPPOV, V.L. POPOV 

fact that the statistics of earthquake force and their time correlations meet the laws of 

Gutenberg-Richter [1,2] and Omori [1,3], typical for self-organized critical systems [4,5], 

is often the basis for the conclusion that earthquakes are a process formed on all spatial 

scales ranging from microscopic to continental plate scale. During last decades [6,10] on 

the basis of modeling by movable cellular automatons and full-scale experiments it was 

shown some possibility of initiation of the fault wing displacement and, consequently, the 

release of a part of the accumulated elastic energy as a result of local effects on the fault 

area, including vibration load and watering [11-16]. 

One of the often applied to study the problems under interest is the well-known 

Burridge-Knopoff (BK) model [17,18] initially proposed many years ago to investigate 

statistical properties of earthquakes. Numerical studies by Carlson et al. [19,20] have 

demonstrated that BK model can reproduce characteristic empirical features of tectonic 

processes such as the Gutenberg-Richter law for the magnitude distribution of 

earthquakes, or the Omori law for statistics of aftershocks [19-22], both properties 

stemming from the so called “self-organized criticality” of this system. It has been 

intensively used to simulate different aspects of the problem [21-40] and to discuss 

general properties of earthquakes statistics as well as predictability of earthquakes. 

Numerical simulations give evidences that the self-organized criticality and the 

corresponding fractal attractor of the system is closely related to dynamic structures with 

“traveling waves” [23], their ordering and specific “phase transitions” [24] controlled by a 

number of parameters (external driving velocity, springs stiffness, number of blocks, their 

mutual interaction and so on).  It is in particular the dependence of the dynamic properties 

of the BK model on spring stiffness which makes it necessary to introduce the 

generalization of the friction law proposed in the works [41,17]. 

The physical reason for the stick-slip instability in the Burridge-Knopoff model is the 

assumed decrease of the friction force with the sliding velocity [19,20]. Motion of a single 

block with this friction law is always unstable which does not correspond to properties of 

real tribological systems.  The real law of rock friction is more complicated [29,34]. In 

this article we base on the realistic friction laws described in the handbooks [34,35]. The 

main qualitative pictures of a realistic law are: (a) approximately logarithmic increase of 

the static friction force as a function of contact time – the property, found already by 

Coulomb, and (b) a logarithmic dependence of the sliding friction on the sliding velocity. 

Both properties can be described in the frame of “state dependent” friction laws by 

introducing additional internal variables describing the state of the contact [17,41-43] as 

well as shear band propagation effect on the deformation and fracture [44,45], generation 

and propagation of slow deformation in geomedia [46-48].  

Modified version of the BK model with a state dependent friction force, reproduces in 

the simplest way both the velocity weakening friction and the increase of static friction 

with time when the block is not moving. 

2. METHOD  

Let us briefly remind the original BK model. The blocks of mass m are attached to a 

moving surface by springs with stiffness k1 and are coupled to each other by springs with 

stiffness k2. The moving surface has velocity v, and the blocks are in contact with a rough 



 Block-media Motion Lattice and Movable Automata Numerical Models 3 

substrate. The friction force F between the blocks and the rough surface is assumed to 

depend only on velocity. The equations of the BK model can be written in the following 

form: 

    
2

2 1 1 12
2 ( )

j

j j j j j

u
m k u u u k vt u F v

t
 


     


, (1) 

where v=const is the external driving velocity and vj≡∂uj/∂t is an array of individual block 

velocities (j=1, … N). The sliding friction force is supposed to decrease monotonously 

from a constant initial value F0. It is further supposed that the static friction F(vj→0) can 

possess any necessary negative value to prevent back sliding 

 
0

0
, 1 ; / 0

1 2 / (1 )( )

( , 0] / 0

j

jj

j

F
F u t

vF v

u t


       

    


  . (2) 

The parameter α defines a rate of friction decrease when block starts to slide, and σ is the 

acceleration of a block at instant when slipping begins [19,20]. The friction equation (2) 

is an oversimplification of real properties of static and kinetic friction which is now 

explained in detail in the books [34,35]. Dieterich has proposed friction equations with 

internal variables [29,30], where the friction force is dependent upon an additional 

variable that describes the state of the contact zone. This variable is in a sense a “melting 

parameter.” Friction force at non-zero velocity drops down from its initial value due to a 

“shear melting” effect which may have different physical nature [17-20]. When the 

motion stops, the surfaces start to form new bonds and the static friction increases with 

time. 

The dynamics of systems with state dependent friction has been investigated in a 

number of papers [29-37]. Majority of these studies were devoted to the simple one-

particle version of the model, but it was also done for many-body BK model with a state 

dependent friction [17], where the friction force is defined by the additional kinetic 

equation: 

  1 0 2
[ ( )]

j j

j j

F v t
F F v

t


  


  . (3) 

with β2<0, at vj>0, and F(vj)=-∞ in opposite case vj≤0. The parameters β1 and β2 have 

following meaning. When a block starts to slide with vj>0 its friction force monotonously 

decreases from an initial value F0. General time scale of this process is determined by the 

first parameter β1 and an effectiveness of the melting is given by a relation between 

absolute value of negative parameter β2<0 and β1.  It should be noted that the kinetic 

equations for state dependent friction at sliding are not unique. Different forms of this 

equation have been proposed [17,18] leading to qualitatively similar behavior of the 

dynamic friction at vj>0. The above form is just simplest and even linear one. 

The equations of motion in Eqs. (1-3) are nonlinear wave equations. Last few decades 

different variants of such equations have been widely studied, starting from the very early 

implementation of the numerical simulations [38]. In particular, it has been shown that N 

interacting segments of the nonlinear chain form a collective attractor with energy transfer 



4 A.E. FILIPPOV, V.L. POPOV 

performed by nonlinear excitations [39,40]. Depending on the relations between total 

energy and strength of the interaction between the blocks k2, the chain can form nearly 

uniform or strongly non-uniform structures and phase patterns.  

Analogous behavior should be expected in the modified Burridge-Knopoff model as 

well. Instantaneously moving blocks must be involved with stationary ones in the overall 

motion. This causes a “detachment” wave, propagating along the chain. Such waves were 

found and studied in the original BK model [18], and the analogous process exists in the 

modified one. There are two possible types of traveling waves: (a) after a transition time 

all the blocks are provoked to move simultaneously, (b) in steady state some of the blocks 

can be found instantaneously motionless while others move. One would expect the 

reaction of the chain to depend on the strength of the interaction between the blocks. This 

question has been studied for the original BK model in [20], and a specific “phase 

transition” between correlated and uncorrelated behavior has been found.   

Real contact of two surfaces is two-dimensional. Let us generalize the MBK to the 

case. The generalized model is very similar to the one-dimensional model but 

incorporates a 2D array of blocks connected by elastic springs in both directions. All 

other components of the MBK remain unchanged. 

The system of equations of motion takes the form: 

    
2

, ,

2 1, 1, , 1 , 1 , 1 , ,2
4 ( );

j n j n

j n j n j n j n j n j n j n

u u
m k u u u u u k vt u F v

t t
   

 
        

 
  (4) 

  
, ,

1 0 , 2 , ,

, ,

( ( ))
, 0, 0

( )                                       0         

j n j n

j n j n j n

j n j n

F v t
F F v v

t

F v v

  


    


  

. (5) 

where j=1, …, Nx and n=1, …, Ny. Nx and Ny are the numbers of elements in the x- and y-

directions. For general applicability to tribological problems Eq.(4) contains also a 

viscous term η∂uj/∂t. Our calculations show that attractor properties are weakly influenced 

by the viscous term at sufficiently high velocities vj. For generality, we keep small 

nonzero value η=0.05. Other parameters are: β1=1, β2=125, F0=1, k1=1. k2=36,  v=0.2. 

Large scale structure of the wave state in the two-dimensional model with number of 

dimensionless grid cells (Nx = 700, Ny = 300) shown by the velocities v=v(x,y;t)≡∂u/∂t is 

reproduced in Fig. 1 by the colored surface (MatLab color-map ‘jet’). Dark blue 

corresponds to zero. Mutual scattering manifests itself in high sharp peaks of the “events” 

reproduced here by the intensive red of the corresponding densities. Instant structure of 

the wave state and history of the events are combined in Fig. 2. The subplots (a) and (c) 

here represent snapshots of the instantaneous densities of the local displacements 

u=u(x,y;t) and velocities v=v(x,y;t)≡∂u/∂t, respectively. Time-space representation of this 

process is shown in the subplot (b), and time dependency of total area of events (d). The 

information about attractor behavior is accumulated in Fig. 3 and briefly described in the 

caption to this figure. 



 Block-media Motion Lattice and Movable Automata Numerical Models 5 

 

Fig. 1 Large scale structure of the wave state in the two-dimensional model shown by the 

velocities v=v(x,y;t)≡∂u/∂t reproduced by the colored surface (MatLab color-map ‘jet’). 

Dark blue corresponds to zero. Mutual scattering manifests itself in high sharp peaks of 

the “events” reproduced here by the intensive red of the corresponding densities. 

 

Fig. 2 Instant structure of the wave state and history of the events in the model. Subplots 

(a) and (c) represent snapshots of the instantaneous densities of the local displacements 

u=u(x,y;t) and velocities v=v(x,y;t)≡∂u/∂t, respectively. The colormaps are equidistant and 

normalized to the intervals between maximum (red) and minimum (blue) values for the 

displacements [-1,1] and velocities [0 1.5]. Time-space representation of this process is 

shown in the (b) and (d) by the spatial-temporal map and time dependency of total area of 

events.  



6 A.E. FILIPPOV, V.L. POPOV 

 

Fig. 3 The attractor of the variable effective friction in projection to the plane {v, Ffriction}, 

where <Ffriction >=<F[vj(t)]+ ηvj(t)>  is shown by color-map ‘jet’ in the subplot (a). The 

attractor behavior of the system causes quasi-periodic behavior of the mean area of the 

events. Correlation function in the subplot (b) here is calculated for the area of events 

shown in Fig. 2a. The same attractor is projected onto the plane {v,u} in subplot (c). 

To characterize the difference of the correlated and uncorrelated behavior, an order 

parameter has been introduced in ref. [25]. If we define a value hj in block number j by 

the condition: 

 
1 u / 0 (the block is moving);

0 u / 0

j

j

j

t
h

t

  
 

  

, (6) 

then the local density of the order parameter H*j can be written as H
*
j=hj(hj+1+hj-1). This 

function takes on unit value H*j=hj(hj+1+hj-1)=1 if block j is moving and exactly one of its 

nearest neighbors is also moving. Further, H*j=hj(hj+1+hj-1)=2 when both nearest 

neighbors of the moving block are in motion. All other cases yield H*j=0. Our 

observations with the modified model show that even for blocks at vanishing inter-block 

interaction k2→0 (when the motion of neighboring blocks is almost uncorrelated), the 

fraction of configurations where sets of neighboring blocks are moving is still relatively 

high. All that to say, the combination H*j=hj(hj+1+hj-1) does not vanish in such an 

uncorrelated system. It is therefore convenient to construct another simple combination: 

 
1 1

1 all 3 blocks  , 1, 1 in contact are moving

0 otherwise
j j j j

j j j
H h h h

 

 
  



. (7) 

This makes Hj=hjhj+1hj-1 zero in all cases except when both neighboring blocks of a 

central sliding block are also in motion. This combination can be used as an “order 

parameter”. To extract integral quantitative information about the steady ordered and 

disordered states we calculate also the time evolution of the ensemble-averaged order 

parameter:  H(t>≡ <Hj(t)>=1/t ∫Hj(t)dt. 

As an illustration the dependence of the integrated order parameter on the stiffness k2, 

showing a transition from correlated to uncorrelated block motion was calculated for one-

dimensional system. It is shown in Fig.4. Two well defined limiting cases can be seen 

easily: the order parameter tends to a constant non-zero asymptote <Hj(t)>→const≠0 at 

strong interaction k2>>1 and vanishes <Hj(t)>→0 below the transition point k2≈0. 



 Block-media Motion Lattice and Movable Automata Numerical Models 7 

 

Fig. 4 Phase transition from correlated to uncorrelated motion of blocks in a one-

dimensional system. The order parameter tends to a constant asymptote at high mutual 

interaction k2>>1 and vanishes below transition point k2= k2,critical≈1. Two intermediate 

fluctuation regions A and B correspond to the states with short-range and long-range 

correlated nonlinear excitations, respectively. 

Let us return now to general problem. If the interaction is already varied due to the 

motion with nonzero velocity, it sounds natural to modify the elasticity (potential channel 

of the interaction) of the system into the state dependent one, which varies at nonzero 

velocity and gradually returns to the initial one in static limit: 

    
2; ,

2
2; 1, 2; 1, 2; , 1 2; , 1 2; , 1 20 2 ,

4
j nj n j n j n j n j n k j n

k
k k k k k k k v

t
   


       


   . (8) 

The same note can be done about the thermalization of the system by a heating produced 

by the mechanical work performed by the external force. Local power consumed by the 

system in every node of the lattice can be calculated as a work per unit of time:  

Q(rj,n)=fj,nvj,n. It can be treated as local source of energy acting in every node and 

temperature expansion inside the system should be described by the thermal conductivity 

equation: 

  , 1, 1, , 1 , 1 , ,4 ( )
j n

j n j n j n j n j n j n

T
T T T T T Q r

t
   


     


 . (9) 

It is possible to repeat all the simulations of the previous section, reproduce all the 

results presented in Figs.1-3 and demonstrate that these properties are common for all the 

modifications of the model. So, let us concentrate below on some additions coming from 

generalized version of the model. 

Large scale structure of the velocity fluctuations analogous to that plotted in Fig. 1 

accompanied by the variations of the elasticity self-consistently calculated according to 

Eq. (8) is shown by the colored surfaces (mainly red picks and blue smooth relief, 

respectively) in Fig. 5. One can see directly that strong fluctuations of the velocity cause 

deep valleys in the surface of elasticity and in turn “prefer” to appear in the same places 

again. 



8 A.E. FILIPPOV, V.L. POPOV 

 

Fig. 5 Large scale structure of the velocity fluctuations accompanied self-consistently by 

the variations of the elasticity are shown by the colored surfaces (mainly red picks and 

blue smooth relief, respectively). 

 

Fig. 6 The temperature deviations (bright spots in Matlab color-map ‘hot’ normalized in 

the interval [0 2]) caused by the same instant configuration of the system as in Fig. 5. 

Green contour lines show median value of the elasticity profile and mark mutual 

correlation between its deep valleys and hot regions. 

The same correlation exists between the spatially dependent elasticity and the 

temperature deviations. In Fig. 6 the temperature deviations (bright spots in Matlab color-

map ‘hot’) caused by the same instant configuration of the system as in Fig. 5 are shown. 

Green contour lines show median value of the elasticity profile and mark mutual 

correlation between its deep valleys and hot regions. 

To extract physically meaningful quantitative information about the earthquakes in the 

system different time moments one can calculate a density of power consumed by the 



 Block-media Motion Lattice and Movable Automata Numerical Models 9 

system at each time moment and overlap it other densities of the system like: distributed 

rigidity, friction constant and temperature.  The local power in the frames of the model 

can be calculated as a work produced by the forces acting on each lattice node of the 

system pij=fijvij during a unit of time. Calculating spatially distributed power pi=fivi 

enables one to evaluate the degree of localization of excitations and hence the fraction of 

the regions in which the fragments of the lattice move practically as a whole. As at a 

particular site the energy is absorbed and released at different points in time, it is useful to 

characterize these sites by the absolute value of the power, |pi|=|fivi|. However, if formally 

accumulate the statistics of such “events” using all the nodes as independent realizations, 

strongly overestimated value of the small “events” will be obtained. Physically only the 

quakes with a magnitude higher than a critical one would be treated as the observable 

ones. Moreover, if the amplitude of the power exceeds the critical value for a number of 

mutually connected nodes, all of them will be treated by the potential observers as one 

common event with an average amplitude in this connected area. So, we organized the 

procedure in this manner: all the nodes were found with the amplitudes higher than the 

critical one; defined all the connected domains of them and calculated mean amplitudes 

for every such a region. To control visually the above procedure, we plotted these 

domains by the (bright) ‘jet’ colors over the (blue) map of the rigidity for every time 

moment. Typical instant realization of such a picture is presented in Fig. 7. The histogram 

calculated for these events, as it is shown in Fig. 8 in double-logarithmic scale. 

 

Fig. 7 Instant connected regions of the “earthquakes” (bright spots) and valleys of the 

small rigidity (deep blue) plotted together to illustrate their mutual correlations. The 

artificial colors for the magnitude and rigidity are shifted to red and blue parts of the 

spectrum to clearly overlap both densities in one picture.   

As we already noted above the variations of the friction, elasticity and local velocities 

are mutually correlated. In other words, new events are preferably generated in the same 

regions (as it is in the reality). As a result, the system partially memorizes such a common 

relief for extremely long times. One can calculate, for example, time-averaged relief of 

effective rigidity, and we did this. It is quite expected that the variations are found to be 

much smaller than local instant variations of the same value, but they exist. 



10 A.E. FILIPPOV, V.L. POPOV 

 

Fig. 8 Scaling of the probability to find an earthquake with a particular value of 

magnitude calculated for the connected event regions. Physically it corresponds to a 

selection of the rare but intensive “events”, which is compatible with the ideology of the 

empirical Gutenberg-Richter law. 

There are also some correlations between such a memorized structure and preferable 

channels of the wave propagation. Again, the correlation is much weaker than it is 

observed between instant variations of different densities, but it still exists. These results 

are illustrated in Fig. 9, where instant and almost stationary time averaged elastic 

coefficient integrated over long time periods ((a) and (b) sub-plots, respectively) shown 

by the ‘jet’ color-map. Black contour line depicts mean value of the long-time living 

rigidity surface.  It marks only partial correlations between the averaged and instant 

rigidities, Analogous (only partial) correlation is found between travelling regions of 

maximal power production (deep red spots in subplot (b)) and time-averaged relief of the 

rigidity. 

 

Fig. 9 Instant and almost stationary time averaged elastic coefficient integrated over long 

time periods ((a) and (b) subplots, respectively) shown by the ‘jet’ color-map. Contour 

line marks only partial correlations between the averaged and instant rigidities, Analogous 

(only partial) correlation is found between travelling regions of maximal power 

production (deep red spots in subplot (b)) and time-averaged relief of the rigidity. 



 Block-media Motion Lattice and Movable Automata Numerical Models 11 

3. MOVABLE AUTOMATA MODEL 

Generally, the variant of the model which uses movable system is based on a 

description of a system of interacting N “particles”. In contrast to described above lattice 

model which treats system as deformable elastic medium with varied friction and rigidity, 

here each “particle” simulates a complete block (“solid plate”) in some instant position 

represented radius-vector ri of its center and the mechanical momentum pi. In simplest 

approximation, the plates interact with a potential U(ri-rj) depending on the distance |ri-rj| 

between their centers. It formally leads to the following many-body Hamiltonian: 

 2

1 , 1

( , ) / 2 (| |) / 2
N N

i i i i i j

i i j

H m U
 

   r p p r r . (10) 

In a static configuration the plates are supposed to occupy some equilibrium positions 

defined by a competition of the repulsion and attraction forces acting between them. In 

the context of this study, it is convenient to represent the interaction potential by a pair of 

the Gauss potentials: 

 2 2(| |) exp{ [( ) / ] } exp{ [( ) / ] }
i j ij i j ij ij i j ij

U C c D d      r r r r r r , (11) 

where Cij and Dij define the magnitude, while cij and dij the radii of attraction and 

repulsion, respectively. The minimization of the energy corresponding to a desired 

equilibrium with some characteristic distance between the centers corresponds to the 

inequalities between the parameters: Cij>>Dij, cij<dij. Initially, correct equilibrium 

positions of the plates are not known and should be found by a self-assembling transient 

process. To obtain it in numerical experiments, the “particles” are initially placed at 

random on a rectangle with the side lengths of [0, Lx] and [0, Ly]. During the transient 

process the system finds some frozen structure which is a compromise between the 

interactions and boundary conditions. 

Depending on the particular study, the boundary conditions can be chosen differently.   

For example, for the case of shear deformation, it is convenient to employ periodic 

boundary conditions on the horizontal axis. It allows simulate potentially infinite run for 

relatively small limited numerical arrays. According to the periodic boundary conditions, 

a particle leaving the interval [0, Lx] is returned to it at the opposite end and the vectors ri-

rj connecting particles either located within the interval [0, Lx] or the images of “escaped” 

particles at the opposite side of the system. On the vertical axis such a system can be 

limited by the horizontal plates at y=0 and y=Ly with reflecting boundary conditions: 

Uup=C·exp[(y-Ly)/c] and Udown=C·exp[-y/c]. Such a combination of the boundary 

conditions simulates together quasi-infinite motion along horizontal axis and confined 

motion between two “external” plates in orthogonal direction. The stronger the 

inequalities C>>Cij and c<<cij, the more rigid and sharp are the walls in relation to other 

forces and lengths relevant to the system. 

Keeping in mind typical “tectonic” geometry with the moving continents, in the 

present study we are interested in a configuration when along horizontal axis the system is 

also limited and moreover, is pressed by a movable “rigid wall”. In analogy with standard 

Burridge-Knopoff model, the plate is driven by an external “elastic” force: 



12 A.E. FILIPPOV, V.L. POPOV 

 
0

( )
ext x

F K L V t  , (12) 

In this case it is natural, in contrast to the periodic one, to take the system which is limited 

along horizontal axis with reflecting boundary conditions: Uright=C·exp[(x-Lx)/c] and 

Uleft=C·exp[-x/c], but with the movable right boundary Lx= Lx(t). 

Since the parameters Cij, Dij, cij and dij form arrays covering all possible combinations, 

it can be said without loss of generality that the Hamiltonian given by Eqs. (10)- (12) 

describes systems with any number of components (plates of the different sizes). What is 

required is just to specify the parameters Cjk and cjk. In a real system different values of all 

the parameters (and what is even more important, the mutual relations between them) 

regulate different sizes of the physical plates and their mutual rigidity. Below, to conserve 

a generality and get statistically representative results, we use the values Cjk, and cjk 

randomly varied with 10 times difference between maximal and minimal ones around unit 

Cjk, cjk=1, and reset them from one run to another of every numerical experiment 

performed at other fixed or systematically varied parameters. 

The simplicity of the equations of motion 

 / ( , ) /
i i i i i i

m t H     v x p p f , (13) 

is deceptive. In majority of the cases, simulation using such equations is extremely time-

consuming because it involves a summation at every time step over all possible neighbors, 

including very distant ones (and even fictitious image particles in the case of periodic 

boundary conditions). However, in the case of the solid plates contacting by their 

boundaries the number of the mutual interactions can be essentially reduced. The 

reduction is based on the following considerations. 

It is natural for this particular problem to construct Voronoi diagram of the set of 

effective “particles” and associate each cell of this diagram with some of the mutually 

contacting plates.  According to mathematical definition, a Voronoi diagram is a partition 

of a plane into regions close to each of a given set of objects. In the simplest case, these 

objects are just set of points in the plane (called seeds). For each seed there is a 

corresponding region, called a Voronoi cell, consisting of all points of the plane closer to 

that seed than to any other. In our case the central “particles” of every “plate” play a role 

of such mathematical “seeds”. 

So, the cells of Voronoi diagram give practically the same picture, which we have in 

mind when building our dynamic model. Every such a cell is surrounded by the limited 

(and in fact quite small) number of the other plates (nearest neighbors). Mathematically 

such neighbors are completely defined by the Delaunay triangulation, which is dual to the 

Voronoi diagram. According to its definition, Delaunay triangulation for a given set P of 

discrete points in a general position is a triangulation DT(P) such that no point in P is 

inside the circumcircle of any triangle in DT(P). As result one gets quite realistic 

representation of the movable plates and short list of the neighbors interacting with every 

particular plate. 

Interacting plates exchange momentum, and hence a dissipation channel acting to 

equilibrate the relative velocities of particles (that happen to be within the distance cv 

close to the equilibrium one) needs to be introduced. This can be done by including an 

additive velocity-depending force. 



 Block-media Motion Lattice and Movable Automata Numerical Models 13 

 2

1

( ) exp [( ) / ]
neigboursN

i i j i j v

j

f c


     
v

v v r r , (14) 

acting on the i-th particle, with some dissipation constant η. As a result, the equations of 

motion assume the following form: 

 i
i i i

m
t


 



r vv
f f . (15) 

They can be integrated by using Verlet’s method, which conserves the energy of the 

system at each time step, provided no energy is supplied externally through mechanical 

work of the movable right boundary which is driven by a balance of the external force 

Fext=K(Lx-V0t) and total interaction with the internal plates repulsed by the boundary 

potential Uright=C·exp[(x-Lx)/c]   

 
2

2
1

( ( )) /
N

x
ext right j x j

j

L
M F U x L t x

t 


    


 . (16) 

So, the equations (15) and (16) must be solved self-consistently. Moreover, counting 

on the nature of the problem under consideration, this system has to be accompanied by 

the equations regulating variation of the effective mutual friction between the plates   with 

negative derivative in the same spirit as it was done above for the standard BK model on 

the elastic lattice: 

  1 0 2( , ) ;
j

j j j

F
F F v t v

t


  


   2 ,<0,  0j nv   and ( )   at 0j jF v v   . (17) 

Taking into account an experience which obtained in modification lattice BK model it 

seems natural to add here also variable intensity of the interactions between the plates at 

nonzero mutual velocities: 

  1 1 10 1 2( , ) ;
j

j j j

C
C C v t v

t


  


   2 ,<0,  0j nv   and 1 10( )   at 0j jC v C v  . (18) 

At the same time, here is no need to include directly the temperature conductivity 

equation into the model because it is already based on the Brownian dynamic equations 

for the plate centres by itself, and in fact automatically involves an exchange by kinetic 

energy (temperature) between them. So, plotting the kinetic energy of the “particles” by 

the artificial colors one gets a color-map in the representation which is natural for the 

model based on the movable automata. The results obtained in the frames of this model 

are summarized in the Figs. 10-15 below. 

Typical instant configuration of the plates pressed from one of the boundaries by an 

external force is shown in Fig. 10. Let us remind that the rigid wall presses the system 

from the right side. As a result, it becomes more compressed form this side and areas of 

the Voronoi cells here become smaller than in average inside the system far from this 

boundary. The mass of every plate conserves in the frame of the model. The only are of it 

changes. Gray-scale map of the subplot (a) is used to show the individual areas of the 

plates.  



14 A.E. FILIPPOV, V.L. POPOV 

The compression and motion of the plates causes nonzero mutual velocities of every 

plate in relation to their proximities. This fact is well seen in the subplot (b) by the 

MatLab ‘jet’ color-map. According to Eq.(17) nonzero mutual velocities result in a 

variation of the effective friction between every plate and its surrounding. This 

characteristic of the system is depicted in the subplot (c). It is important to note also that 

the Delaunay triangulation is recalculated at every time step. As a result, the links 

between the plates moving with different velocities can appear and disappear depending 

on their physical positions at a particular time moment. However, in sufficiently dense 

system (which is used here) such changes of the neighbors happen quite rarely. So, the 

model seems quite realistic to describe some features of the tectonic motion.   

 

Fig. 10 Typical instant configuration of the plates pressed from one of the boundaries by 

an external force. Gray-scale map of the subplot (a) shows the areas of the plates more 

compressed from one side (white color) than inside the system (different intensity of gray 

color). The compression and motion of the plates causes different (mean) mutual velocity 

of every plate in relation to its proximity which results in the variation of the effective 

friction between this plate and surrounding ones. Both these characteristics are plotted by 

the ‘jet’ colors. 



 Block-media Motion Lattice and Movable Automata Numerical Models 15 

Fine structure of the model is presented in Fig. 11. Movable Voronoi diagram in 

subplot (a) here show the plates in the frame of the model, as above. However, here they 

are accompanied by the colored seeds of the Voronoi cells. The colors show different 

instant velocities of the plates. Dual to the above diagram Delaunay triangulation shows 

connections between the plates and their neighbors counted in the model. According to its 

mathematical definition it naturally gives a list of the neighbors for every plate. So, one 

can calculate the interactions between them, solve the dynamic equations, so on. Different 

sizes of the nodes here reproduce different (randomly distributed) interaction parameters. 

Besides, it is convenient to control instant number of them for every cell (symmetry of the 

plate polygon) and visualize it by the colors from the same standard ‘jet’ map. 

 

Fig. 11 Movable Voronoi diagram in subplot (a) describes the plates in the frame of the 

model. The colors of the seeds show different instant velocities of the plates. Dual to the 

above diagram Delaunay triangulation shows connections between the plates and their 

neighbors counted in the model. Different sizes of the nodes here reproduce different 

interaction coefficients and their colors show different instant number of neighbors. 



16 A.E. FILIPPOV, V.L. POPOV 

Two projections of the phase portrait on different planes (subplots (a) and (b)) are 

presented in the Fig. 13 by the colored circles of different sizes. The colors and sizes of 

the circles have the same meaning as in the Fig. 12a.  Green line with the arrows shows 

typical individual trajectory, which every point passes through the phase space with the 

time and reproduces the evolution of the effective friction due to periods of large or small 

local velocities.  Generally, the subplot (b) here has the same physical meaning as the 

phase portrait plotted in Fig. 3a in the lattice variant of BK model. 

 

Fig. 12 The same instant configuration as in Fig. 11 reproducing mutual relations between 

the interactions, local rigidities, velocities and friction coefficient directly marked in the 

plots. The colors and sizes of the circles in the subplot (a) depict velocities and interaction 

coefficients of the individual plates, respectively. Yellow, cyan, magenta and red lines in 

the subplot (a) display integrated over y-coordinate averaged x-dependencies of the 

velocity, rigidity, effective friction coefficient density and temperature respectively. The 

colors of the circles in the subplot (b) correspond to different individual power of the 

events in every plate. 



 Block-media Motion Lattice and Movable Automata Numerical Models 17 

To get the local power in the frames of the discrete variant of the model we calculate 

the work produced by the forces acting on each plate (seed) of the system pi=fivi. 

Calculation of spatially distributed power pi=fivi enables one to evaluate the degree of 

localization of excitations and hence the fraction of the regions in which the fragments of 

the lattice move practically as a whole. As at a particular site the energy is absorbed and 

released at different points in time, it is useful to characterize these sites by the absolute 

value of the power, |pi|=|fivi|. It is shown in subplot (b) of the Fig.12. 

Summation over all particles yields the total power VFext(t)=P(t)=∑fivi. It gives an 

energy which is consumed during every time unit due to motion of rigid right boundary, 

which is forced by external flow with constant velocity, V=const.  However, the boundary 

here is moving with some varied velocity V(t) which is determined by a balance of the 

total force, Fext(t) and only its time average coincides with <V(t)>=V0 . 

 

Fig. 13 Two projections of the phase portrait on different planes (subplots (a) and (b)). 

The colors and sizes of the circles depict velocities and interaction coefficients of the 

individual plates, respectively. Green line with the arrows shows typical individual 

trajectory, which every point passes in phase space with the time. 



18 A.E. FILIPPOV, V.L. POPOV 

 

Fig. 14 Instant configuration of the density (surface in subplot (a)) and time-space 

diagram of the history of events (colormap in subplot (b)). Bold black contours mark 

intensive “earthquake” events near the propagating right boundary. Different speed of the 

front and some “epochs” of intensive events can be seen from different inclination of the 

front and horizontal lines in the recorded history. 

In particular, it causes quasi-periodic oscillations of the velocity of boundary which 

are typical for the “stick-slip” nature of the processes. It is actually the same phenomenon, 

which takes place in ordinary BK model, shown above in Figs. 2-3. In the discrete model 

one can additionally observe the waves of the velocities, densities and power of the 

events, which are accompanied by quasi-periodic “epochs” of intensive “earthquakes”. 

Instant configuration of the density and history of the density integrated along y-axis 

corresponding to the magenta curve in Fig. 12a are shown in Fig. 14. The instant surface 

is plotted in subplot (a) and corresponding time-space diagram of the history of events is 

presented in the subplot (b). Bold black contours mark intensive “earthquake” events near 

the propagating right boundary. Different speed of the front and some “epochs” of 

intensive events can be seen from different inclination of the front and deeply colored 

almost horizontal lines in the recorded history. 

Exactly the same history of the events (recorded during the same calculation run) as in 

Fig.  14, is depicted by the intensity of power production (shown by the ‘hot’ color map). 

The epochs of the intensive earthquakes are much stronger pronounced here than in the 

previous figure. 



 Block-media Motion Lattice and Movable Automata Numerical Models 19 

 

Fig. 15 Exactly the same history of the events, as in Fig.14, depicted by the intensity of 

power production (shown by the ‘hot’ colormap). The epochs of the intensive earthquakes 

are much stronger pronounced here than in the previous figure. 

4. CONCLUSION 

It is shown that original Burridge-Knopoff model allows a number of the 

modifications which can increase its ability to be adapted to description of the different 

practically important features of the tectonic phenomena. Basic modifications of the 

model were already made previously. They incorporate naturally introduced state 

dependent friction where the friction force is defined by the additional kinetic equation as 

well as extension of the model to more realistic two-dimensional case.  

In the present paper we continued previous modifications and mainly concentrated on 

the following extensions of the model. First of all, elastic force and thermal conductivity 

were added into consideration. Besides to the velocity weakening of friction between 

moving blocks, we combined an increase of both static friction and rigidity during stick 

periods as well their simultaneous weakening during motion. It is expected to give better 

and qualitatively correct description of plural transitions from smooth sliding to stick-slip 

behavior, which are normally observed in various tribological and tectonic systems. To 



20 A.E. FILIPPOV, V.L. POPOV 

get easier comparison with the previous versions of the model we initially did this using in 

standard approach based on interacting elastically connected mesh elements.   

An alternative version of the model is finally proposed.  It applies a simulation of the 

motion of relatively large solid blocks, utilizing technique of the movable cellular 

automata. In this advanced version we also involve all the components of the real system: 

state depending friction, changeable rigidity, as well as heat production and thermal 

conductivity due to dynamics of the interacting plates, which were already included into 

the previous versions based on traditional lattice model.  

Proposed modification has its own advantages and disadvantages, which were 

discussed in the paper step by step. However, a combination of the different approaches 

gives additional degrees of freedom, which allows select more convenient or adequate 

approach depending on the particular problem under consideration. We plan to continue 

further and combine our studies in both alternative directions. 

REFERENCES  

1. Gutenberg, B., Richter, C.F., 1994, Frequency of earthquakes in California, Bulletin of the 

Seismological Society of America, 34(4), pp. 185-188. 

2. Bak, P., Christensen, K., Danon, L., Scanlon, T., 2002, Unified scaling law for earthquakes, Physical 

Review Letters, 88, 178501.  

3. Žalohar, J., Omori's Law, 2018, Chapter 10 in the book “Developments in Structural Geology and 

Tectonics”, Elsevier, pp. 123-134. 

4. Bak, P., Tang, C., Wiesenfald, K., 1987, Self-organized criticality: An explanation of the 1/f noise, 

Physical Review Letters, 59, 381. 

5. Carlson, J.M., Langer, J.S., 1989, Properties of earthquakes generated by fault dynamics, Physical 

Review Letters, 62, 2632. 

6. Ruzhich, V.V., Smekalin, O.P., Shilko, E.V., Psakhie, S.G., 2002, About nature of “slow waves” and 

initia-tion of displacements at fault regions, In: Proceedings of the international conference on new 

challenges in mesomechanics, pp.311-318. 

7. Ruzhich, V.V., Truskov, V.A., Chernykh, E.N., Smekalin, O.P., 1999, Neotectonic movements in fault 

zones of the Baikal region and their origin mechanisms, Russian Geology and Geophysics, 40(3), pp. 

356-368. 

8. Psakhie, S.G., Shilko, E.V., Astafurov, S.V., 2004, Peculiarities of the mechanical response of 

heterogeneous materials with highly deformable interfaces, Technical Physics Letters, 30, pp. 237-239. 

9. Psakhie, S.G., Ruzhich, V.V., Shilko, E.V., Popov, V.L., Dimaki, A.V., Astafurov, S.V., Lopatin, V.V., 

2005, Influence of the state of the inter-block interfaces on the character of local displacements in 

block-structured media, Technical Physics Letters, 31, pp. 712-715. 

10. Ruzhich, V.V., Psakhie, S.G., Bornyakov, S.A., Smekalin, O.P., Shilko, E.V., Chernykh, E.N., 

Chechelnitsky, V.V., Astafurov, S.V., 2002, Investigation of influence of vibroimpulse excitations on 

regime of displacements in seismically active fault regions, Phys. Mesomech., 5(5), pp. 85-96. 

11. Bhattacharya, P., Chakrabarti, B.K., Samanta, K.D., 2009, Fractal models of earthquake dynamics, 

Heinz Georg Schuster (ed), Reviews of Nonlinear Dynamics and Complexity, Wiley-VCH, pp. 107–150. 

12. Ford, J., 1992, The Fermi-Pasta-Ulam problem: Paradox turns discovery, Physics Reports, 213(5), pp. 

271-310. 

13. Lichtenberg, A.J., Lieberman, M.A., 1992, Regular and Chaotic Dynamics, Springer, New York, 692 p. 

14. Fillipov, A.E., Hu, B., Li, B., Zeltser, A., 1998, Energy transport between two attractors connected by a 

Fermi-Pasta-Ulam chains, Journal of Physics A: Mathematical and General, 31, 7719. 

15. Filippov, A.E., Popov, V.L., Psakhie, S.G., 2006, Converting displacement dynamics into creep in block 

media, Technical Physics Letters, 32(6), pp. 545-549. 

16. Filippov, A.E., Popov, V.L., Psakhie, S.G., 2008, Correlated impacts optimizing the transformation of 

block me-dium dynamics into creep regime, Technical Physics Letters, 34(8), pp. 689-692. 



 Block-media Motion Lattice and Movable Automata Numerical Models 21 

17. Filippov, A.E., Popov V.L. 2021, Study of Dynamics of Block-Media in the Framework of Minimalistic 

Numerical Models, Chapter 7 in the book «Multiscale Biomechanics and Tribology of Inorganic and 

Organic Systems», Springer International Publishing, pp. 143-168. 

18. Burridge. R., Knopoff, L., 1967, Model and Theoretical Seismicity, Bulletin of the Seismological 

Society of America, 57(3), pp. 341-371. 

19. Bruce, S.E., Carlson, J.M., Langer, J.S., 1992, Patterns of seismic activity preceding large earthquakes, 

Journal of Geophysical Research, 97, pp. 479-488. 

20. Carlson, J.M., Langer, J.S., Shaw, B.E., 1994, Dynamics of earthquake faults, Review of Modern 

Physics, 66(2), pp. 657-670. 

21. Saito, T., Matsukawa, H., 2007, Size dependence of the Burridge-Knopoff model, Journal of Physics: 

Conference Series, 89, 012016. 

22. Mandelbrot, B., 1982, The Fractal Geometry of Nature, Freeman and Co., San Francisco, 460 p.  

23. Pelletier, J.D., 2000, Spring-block models of seismicity: review and analysis of a structurally 

heterogeneous model coupled to the viscous asthenosphere. Geocomplexity and the Physics of 

Earthquakes, Geophysical Monographs, American Geophysical Union, pp.27-42.  

24. Muratov, C.B., 1999, Trave ling wave solutions in the Burridge-Knopoff model, Physical Review E, 59, 

3847.  

25. Huisman, B.A.H., Fasolino, A., 2005, Transition to strictly solitary motion in the Burridge-Knopoff 

model of multicontact friction, Physical Review E, 72, 016107. 

26. Clancy, I., Corcoran, D., 2005, Criticality in the Burridge-Knopoff model, Physical Review E, 71, 

046124. 

27. Clancy, I., Corcoran, D., 2006, Burridge-Knopoff model: Exploration of dynamic phases, Physical 

Review E, 73, 046115. 

28. Mori, T., Kawamura, H., 2008, Simulation study of earthquakes based on the two-dimensional 

Burridge-Knopoff model with long-range interactions, Physical Review E 77, 051123. 

29. Dieterich, J.H., 1974, Earthquake mechanisms and modeling, Annual Review of Earth and Planetary 

Sciences, 2, pp. 275-301. 

30. Dieterich, J.H., 1978, Time-dependent friction and the mechanics of stick-slip, Pure and applied 

geophysics, 116, pp. 790-806.  

31. Rice, J.R., Gu, J.C., 1983, Earthquake aftereffects and triggered seismic phenomena, Pure and applied 

geophysics, 121, pp. 187-219. 

32. Gu, J.C., Rice, J.R., Ruina, A.L., Tse, S.T., 1984, Slip motion and stability of a single degree of freedom 

elastic system with rate and state dependent friction, Journal of the Mechanics and Physics of Solids, 

32(3), pp. 167-196. 

33. Marone, C., 1998, Laboratory-derived friction laws and their application to seismic faulting, Annual 

Review of Earth and Planetary Sciences, 26, pp. 643-696. 

34. Persson, B.N.J., 2000, Sliding Friction, Physical Properties and Applications, Springer, Berlin. 

35. Popov, V.L., 2010, Contact Mechanics and Friction Physics, Springer, Berlin, 362p. 

36. Popov, V.L., 2000, A theory of the transition from static to kinetic friction in boundary lubrication 

layers, Solid State Communications, 115(7), pp. 369-373. 

37. Filippov, A.E., Klafter, J., Urbakh, M., 2004, Friction through dynamical formation and rupture of 

molecular bonds, Physical Review Letters, 92, 135503. 

38. Fermi, E., Pasta, J., Ulam, S., Tsingou, M., 1993, The Many-body Problem: An Encyclopedia of Exactly 

Solved Models in One Dimension, World Scientific, Singapore, 988 p. 

39. Braun, O.M., Hu, B., Filippov, A.E., Zeltser, A., 1998, Traffic jams and hysteresis in driven one-

dimensional systems, Physical Review E, 58, 1311. 

40. Filippov, A.E., 1994, Mimicry of phase transitions and the large-river effect, JETP Lett., 60(2), pp. 133-

137. 

41. Filippov, A.E., Popov, V.L., 2010, Modified Burridge–Knopoff model with state dependent friction, 

Tribology International, 43, pp. 1392-1399. 

42. Ohmura, A., Kawamura, H., 2007, Rate-and state-dependent friction law and statistical properties of 

earthquakes, Europhysics Letters, 77, 69001. 

43. Persson, B.N.J., Popov, V.L., 2000, On the origin of the transition from slip to stick, Solid State 

Communications, 114(5), pp. 261-266.  

44. Balokhonov, R., Romanova , V., Zemlianov, A., 2021, A mesoscopic analysis of a localized shear band 

propagation effect on the deformation and fracture of coated materials, Reports in Mechanical 

Engineering, 2(1), pp. 6-22.  



22 A.E. FILIPPOV, V.L. POPOV 

45. Romanova, V., Balokhonov, R., Zinovieva, O., 2021, Mesoscale deformation-induced surface 

phenomena in loaded polycrystals, Facta Universitatis-Series Mechanical Engineering, 19(2), pp. 187-

198. 

46. Kocharyan, G.G., Kishkina, S.B., 2021, The Physical Mesomechanics of the Earthquake Source, Phys. 

Mesomech., 24, pp. 343–356. 

47. Makarov, P.V., Khon, Y.A. Autosoliton View of the Seismic Process. 2021, Part 1. Possibility of 

Generation and Propagation of Slow Deformation Autosoliton Disturbances in Geomedia, Phys. 

Mesomech., 24, pp. 363–374. 

48. Makarov, P.V., Smolin, I.Y., Khon, Y.A., et al. 2021, Autosoliton View of the Seismic Process. Part 2. 

Possibility of Generation and Propagation of Slow Deformation Autoso liton Disturbances in 

Geomedia, Phys. Mesomech., 24, pp. 375–390.