On the similarity of shear deformation of a granular massif and a fragmented medium in the seismically active area

In recent research, the dynamics of the medium located in the seismic region at the boundary of tectonic plates is considered as the behavior of a complex open system that is in a state of self-organized criticality. Such an approach results from the very laws of earthquake generation and the complex structure of these areas. The network of faults and cracks makes seismic zones significantly heterogeneous and fragmented. Therefore, discrete models are increasingly used to model the dynamics of these media. The basis for comparing the model and the full-scale object serves the statistical regularities of their dynamic deformation. Relying on this concept, in the paper it is modeled the shear dynamics of a granular massif composed of identical cubic granules and is compared system’s statistical characteristics with the similar characteristics obtained for the earthquake generation zone. Shear deformation is carried out by means of the box consisting of two parts — movable and immovable ones. The movable part possesses the cover which receives kinetic energy from the granular massif in the process of shear deformation. For numerical simulations of the shear dynamics, the discrete element method is applied. The numerical calculations result in the distribution of cover’s kinetic energy jumps simulating the perturbations transmitted from the granular system to an external medium. It turned out that the distribution for these perturbations is the power dependence with an exponent that is inherent in earthquakes (Gutenberg-Richter law). Before and after large perturbations it is observed the swarms of smaller perturbations which are the analogues of foreshocks and aftershocks. The distributions of element’s velocity fluctuations and the correlation of velocity fluctuations are calculated as well. It is revealed the similarity of distributions for velocity fluctuations in the model massif and in the seismically active region of California, which includes the San Andreas fault. Moreover, the similarity of corresponding correlation functions is shown. They both are the functions of the stretched exponent. The obtained result indicates that shear processes in granular massifs and natural seismic processes in the San Andreas Fault are statistically similar.

occurs intermittently and it is accompanied by the generation of acoustic perturbations. These stochastic perturbations obey the same statistical laws as earthquakes: Gutenberg-Richter law, Omori's law, the universal low for distribution of intervent times between avalanches [Geller et al., 2015;Lherminier et al., 2019;Mykulyak et al., 2019bMykulyak et al., , 2021Kumar et al., 2020].
Meroz and Meade studied in detail the motion of tectonic structural elements in California, in the area between two tectonic plates: the Pacific and the North American, which includes the San Andreas Fault [Meroz, Meade, 2017]. This area is divided by a network of faults into a large number of tectonic elements, which as a result of the tectonic plate movement are also involved in the movement. Their studies have shown that the behavior of this area is similar to the shear deformation of the granular medium. In this case, the distribution of velocity fluctuations differs from the Gaussian distribution, revealing «heavy» tails, and the correlation function attenuates as a stretched exponent [Meroz, Meade, 2017].The authors compared the velocity field with the two-dimensional velocity field, which takes place in the process of shear deformation of the granular medium with spherical grains.
In this paper, we study the properties of the velocity field in a three-dimensional model granular medium formed by cubic elements. This problem statement is closer to the real shear process in the seismically active zone in comparison with the problem of shear of smooth spherical granules because in the kinematics of the granular medium a significant role is played by the rotational degrees of freedom. Moreover, the three-dimensional shear process is more in line with the natural process.
The model for simulating the shear dynamics of cubes. The simulation of shear granular dynamics is performed by the discrete element method (DEM) [Cundall, 1971]. The well-known advantage of this method consists in the ability to describe the dynamics of each element of structure in detail.
The equations of motion of the structural elements have the following form: Here r ci is the coordinate of the center of i-th cube and m i is its mass; F ij is the contact force appearing between i-th and j-th cubes; K ci is the kinetic moment of i-th cube with respect of its center; M c (F ij ) is the torque with respect to the center of i-th cube. N is the total number of elements. The summation is performed for all j-th cubes contacting with i-th cube.
The kinetic moment is related to the body's angular velocity via the tensor of moments of inertia. In the general case, their directions do not coincide. Due to the symmetry of the cube, its central ellipsoid of inertia degenerates into a sphere. Therefore, the vectors of the kinetic moment and angular velocity always have the same direction, i. e. K ci =I ci w ci , where w ci is the vector of angular velocity for the cube rotating about the center, The description of the translational motion of the cube center is a relatively simple task, whereas the consideration of rotational motion demands more effort. The use of the quaternion technique [Gürlebeck, Sprössig, 1997] greatly facilitates the task. In the motionless reference frame, the kinematic equation for the rotation quaternion has the fol- λ 2 , λ 3 )=λ 0 +λ 1 i 1 +λ 2 i 2 +λ 3 i 3 =λ 0 +λ is the rotation quaternion. During calculations, we also perform the correction of the quaternion norm to 1 in accordance with the modified equation The parameter q belongs to the interval (0; 1). We put q=0.5. Thus, relations (1)-(2) form the closed system of equations for describing the interaction of pair of cubes. The transformation of the vector r to the vector r′ during rotation, that is described by the quaternion Λ, is determined by the formula ′ = Λ Λ r r   . To solve Eqs. (1), Beeman's algorithm is used [Beeman, 1976]. Equation (2) is solved by the Runge-Kutta 4-th order method.
The interaction between bodies is determined by a single force, which decomposes into a normal force F n and a tangential force The direction of the normal force F n is defined by the contact normal. To derive the contact normal, we use the vertices of the contact line which is the broken line of intersection for the surfaces of two contacting particles [Zhao et al., 2015]. Based on these vertices, we construct the fitting plane by means of the linear square fitting. To do this, the singular value decomposition method [Shakarji, 1998;Forsythe et al., 1977] is applied. The normal of this plane coincides with the contact normal and defines the direction for the normal force of interaction between particles. To determine the absolute value of the normal force, a model proposed by Nassauer and Kuna [Nassauer, Kuna, 2013] is used that is suitable for finding this force for particles of different shapes. According to this model, the contact interaction force between bodies depends on the volume V of the intersection domain of interacted bodies and the depth of penetration d: and d max is the maximal distance between vertices of the overlapping domain (direction l m ), d p is the maximal distance between vertices in the direction n p =n f ×l m . The depth of penetration d is defined as a distance of penetration of intersection domain in the direction of normal force and is derived by the following formula: d=max(n f •p c )-min(n f •p c ), p c is the vertex coordinates of the penetration figure.
The tangential forces F f are determined using a friction model taking into account static friction (adhesion) F s <m s F n and kinetic friction (sliding) F k <m k F n , where m s and m k are the empirical coefficients. The transition from one model to another one occurs smoothly using the ratio proposed in [Nassauer, Kuna, 2013]: Here v t is the tangential velocity, v s is the velocity of transition from static to kinetic friction. The direction of friction force F f is parallel to the relative tangential velocity of cubes at the contact point. The application of this algorithm results in finding all pairs of contacting cubes and the characteristics of their contact. Based on the described algorithm, the code cuBluck has been developed [Mykulyak et al., 2019а].
In the simulations, it is fixed the following constants: density r=1.2•10 3 kg/m 3 ; Young's modulus E=3•10 9 N/m 2 ; Poisson's ratio v=0.3; static friction coefficients m s =0.7; kinetic m k =0.4; the characteristic velocity of transition from static to kinetic friction is v s =0.01 m/s. Granular medium shear at a constant rate of deformation. Shear deformation of the granular massif is carried out due to the movement of the upper part of the box. It consists of two parts, namely, the fixed lower part with a height of 60 mm and the movable upper part of height of 80 mm. In total, the massif contains 3000 cubes of the same size 10 mm. At the top it is located the cover weighing 14 kg, which can move freely in the vertical direction. The bottom has the size 0.2×0.3 m. The box filled with randomly placed cubes is shown in Fig. 1. The massif used for modeling is obtained with compacting cubes randomly placed initially in the volume of 0.2×0.3×0.5 m.
The velocity of the upper box part is constant V b =1 м/с. For clarity, the locations of the cubes at two different moment of time are presented in Fig. 2. In the process of shear deformation, it is occurred dilatation causing the cover's vertical movement. The movement of the cover is intermittent in the form of a sequence of jumps. The dependence of the cover's kinetic energy on time is shown in Fig. 3, a. The fragment of this dependence shown in Fig. 3, b indicates that this curve has many local maxima and minima. Energy jumps from adjacent minima to maxima can be considered as excitations that are transmitted from the granular system to the environment. The total number of perturbations is N e =7464. This number is sufficient to investigate the statistical properties of the perturbation sequence. Fig. 4 shows the cumulative complementary distribution of radiated energy jumps DE c in the logarithmic coordinates. The straight line has a slope b=0.94±0.01.This model reproduces accurately the Gutenberg-Richter law   On the similarity Of shear defOrmatiOn Of a granular massif and ... Геофизический журнал № 3, Т. 43, 2021 165 in the energy representation ( ) ∝ , in which the exponent varies in the range 0.80-1.05 [Olami et al., 1992].
Statistical analysis of the sequence of perturbations shows that it is possible to identify large jumps before and after which the swarms of smaller amplitude perturbations are observed. This is evidenced by the dependence of the average relative number of disturbances before and after the main disturbance shown in Fig. 5.
Fluctuations of element's velocities in the shear process. Meroz and Meade [Meroz, Meade, 2017] studied in detail the velocity field observed over a period of about 10 years in the area between the Pacific and North American plates in California. They examined the area from the Pacific Plate, where they determined y=0, to the North American Plate y=565 km, which runs parallel to the San Andreas Fault. The axis x is directed parallel to the San Andreas Fault. They have calculated then the average velocity profile along the y axis, ( ) i y v , by discretizing the y axis into n=65 layers, with dy=8.7 km. Velocity fluctuations in each layer were also calculated as ( ) is the average velocity in the i-th layer. It turned out that the fluctuation distribution calculated later is significantly different from the Gaussian distribution. To characterize the system's kinematics, which is determined by velocity fluctuations, Meroz and Meade also calcu-lated the correlation function were r is a distance between two reference points. The authors found that the correlation function is best approximated by a stretched exponential function where values b=0.75, and x=92 km obtained by the method of least squares.
In this paper we compare the above statistical characteristics of tectonic movements with similar characteristics obtained in the simulation of 3D shear deformation of the granular array with cubic grains. For this purpose, the entire volume of the granular massif is divided into 20 equal horizontal layers 6 mm thick and the average velocity < v y > is calculated in each layer. The dependence of the average velocity on the z coordinate is shown in Fig. 6. It is seen that the elements at the bottom of the box are motionless, and near the cover move at the velocity of 1 m/s coinciding with the cover velocity. Between these areas there is a transition area where the average velocity varies. For this region, namely for 0.04 < z < 0.10, fluctuations are calculated and the distribution function is constructed (Fig. 7, a). This distribution, as in [Meroz, Meade, 2017], is best approximated by an exponential function Φ(r)=Ae -r/z .  In our case, the coefficients are as follows: A=473±12, z=0.13±0.02. We also have built correlation function (3), which is shown in Fig. 7, b. Approximation of this dependence by the function of the stretched exponent (4) using the method of least squares gives the following values of constants: b=0.96±0.05, x=0.014±0.001.
Thus, the velocity field in shear deformation of the model granular medium has similar properties to the velocity field in the boundary region, namely, the distribution of fluctuations relative to the mean velocity also has an exponential dependence, not Gaussian, and the spatially correlated fluctuation function is best approximated by the stretched exponent higher than the transition zone in the San Andreas fault area.
Concluding Remarks. The simulation results confirm the complexity of the behavior of the granular medium formed from ribbed, namely cubic, grains. This complexity is similar to the complexity that is inherent in the natural media in seismic areas. We obtained that the shear deformation of the granular medium has an intermittent stochastic character. The vertical movement of the cover, which is caused by the dilatancy effect, also occurs intermittently. The distribution of the cover's kinetic energy jumps manifests the power nature, similar to that for earthquakes. Jumps of kinetic energy can be considered as the energy of perturbations that radiates the granular array into the environment. In addition, both aftershocks and foreshocks occur for large energy spikes, like in natural seismic processes.
The distribution of velocity fluctuations is also constructed and the correlation of velocity fluctuations is calculated. The similarity of the distributions of velocity fluctuations in the model medium and in the seismically active region in California, which includes the San Andreas fault, is revealed. There is also a similarity of correlation functions: in both cases they are functions of the stretched exponent. This analogy makes it possible to investigate in more detail the properties of such complex phenomena as earthquakes.
The work is partially supported by the projects 0118U000043 and 0118U000044.