Laboratory Observations of Linkage of Preslip Zones Prior to Stick-Slip Instability

Field and experimental observations showed that preslip undergoes a transition from multiple to single preslip zones, which implies the existence of linkage of preslip zones before the fault instability. However, the observations of the linkage process, which is significant for understanding the mechanism of earthquake preparation, remains to be implemented due to the limitations of observation methods in previous studies. Detailed spatiotemporal evolutions of preslip were observed via a high-speed camera and a digital image correlation method in our experiments. The normalized length of preslip zones shows an increase trend while the normalized number of preslip zones (NN) shows an increase followed by a decrease trend, which indicate that the expansion of the preslip undergoes a transition from increase to linkage of the isolated preslip zones. The peak NN indicates the initiation of the linkage of preslip zones. Both the linkage of the preslip zones and the decrease in the normalized information entropy of fault displacement direction indicate the reduction of spatial complexity of preslip as the instability approaches. Furthermore, the influences of dynamic adjustment of stress along the fault and the interactions between the asperities and preslip on the spatial complexity of preslip were also observed and analyzed.


Introduction
The reduction of the disaster and losses caused by destructive earthquakes is a matter of concern to all humans. One fundamental but difficult solution is improving the accuracy of earthquake prediction based on the observations of earthquake precursors. Studying the physical process of earthquake preparation is the basis for understanding the mechanism of earthquake precursors. The motion of a seismogenic fault is jerky, which can be well described as "stick-slip" [1]. The heterogeneities in material composition, geometrical morphology, and mechanical properties of the fault result in the occurrence of its overall instability preceded by local instability. The local instability, which is caused by preslip in pre-existing faults, establishes a bridge from the stable "stick" state to the unstable "slip" state of the fault. The occurrence of preslip prior to the overall instability of a fault was found prevailing in laboratory [2,3], which is an important physical process during the earthquake preparation of faults. A large number of experiments have shown that the fault experiences an acceleration process in the preslip speed as well as the preslip zone expansion speed [3][4][5][6][7][8][9].
An earthquake nucleation model based on fracture mechanics or friction rate and state dependent constitutive relation [10,11] was proposed and stressed that the instability of a fault occurs when the preslip zone expands to a certain critical scale. The migration of foreshock sequences, which represented the propagation front of the preslip zone, was considered to be the earthquake   recorded by the camera. The images were processed by a DIC method [9,18,[29][30][31][32][33][34] to calculate the deformation of the sample. Two symmetric measurement lines (l1 and l2 in Figure 1a), each of which was composed of 1761 measurement points (pixels) and was 4.2 mm offset from the fault, were used to calculate the fault displacement and strain. Each couple of symmetric measurement points from the two measurement lines made up a fault displacement sensor, which meant that 1761 fault displacement sensors with spacing of 0.2 mm were used to measure the fault displacement. This dense distribution of fault displacement sensors enabled to observe the detailed spatial evolution of the preslip of the fault. The error of the fault displacement measurement was ±5 μm. For the details of the setting parameters of the DIC method, please see the Section 2.2.  Variation of the friction coefficient (μ) and piston displacement along the Y axis (dy) with time (t). Stick-slip events labeled by E1 and E2 were analyzed in this study. Figure 2. Variation of the friction coefficient (µ) and piston displacement along the Y axis (dy) with time (t). Stick-slip events labeled by E1 and E2 were analyzed in this study.
We used a high-speed camera (FASTCAM SA2, Photron, Tokyo, Japan) to capture images of the upper sample surface during Events E1 and E2 with acquisition durations of 2.8 s and 5.5 s, respectively. The sampling rate was 1000 frames per second. The resolution of each image was 2048 × 2048 pixels, and the actual size of each pixel corresponding to the surface of the sample was 149.5 × 149.5 µm 2 . It took 30 min for the computer to download the image data from the cache of the camera, during which the camera couldn't record new images. Therefore, only Events E1 and E2 were recorded by the camera. The images were processed by a DIC method [9,18,[29][30][31][32][33][34] to calculate the deformation of the sample. Two symmetric measurement lines (l 1 and l 2 in Figure 1a), each of which was composed of 1761 measurement points (pixels) and was 4.2 mm offset from the fault, were used to calculate the fault displacement and strain. Each couple of symmetric measurement points from the two measurement lines made up a fault displacement sensor, which meant that 1761 fault displacement sensors with spacing of 0.2 mm were used to measure the fault displacement. This dense distribution of fault displacement sensors enabled to observe the detailed spatial evolution of the preslip of the fault. The error of the fault displacement measurement was ±5 µm. For the details of the setting parameters of the DIC method, please see the Section 2.2.

Parameter Setting of DIC Method
The DIC method, which is a full-field and noncontact-type deformation measurement method based on pattern matching, has been widely used in experimental mechanics [9,18,[32][33][34]. We also used this method to calculate the deformation field of the sample surface with the digital images of the original speckles of the sample surface obtained in the present experiment. The factors that influence the accuracy of the DIC measurement include the hardware system, the environment, the gray-level distribution in the image, the parameters set for calculation, etc. The systematic strain measurement error induced by self-heating of the cameras is remarkable and should be taken into account before measurement. According to a test of six cameras of four different types, a preheating period of at least 1-2 h is required to reduce the systematic strain measurement error [33]. In the present experiment, the preheating time was 3 h.

Influence of the Subregion Size on the Correlation Coefficient
A subregion centered at a point in the original image is chosen to calculate the correlation coefficient (CC) between the subregion and the moving subregion of the same size in the deformed image when the location of the point in the deformed image is searched. CC is an important parameter reflecting the subregion matching quality and the accuracy of the DIC calculation. CC ranges from 0 to 1 and increases as the subregion matching quality improves. It also depends on the subregion size ( Figure 3). The average CC attains a peak value when the radius of a subregion R is 19 pixels. Although the standard deviation of CC decreases as R increases, the decrease of the standard deviation of CC is remarkable when R is less than 19 pixels but is inconspicuous when R is greater than 19 pixels. On the other hand, because the computational time and space filtering effect increase as the size of a subregion increases, a smaller subregion results in a better DIC calculation under the premise of a desirable CC. Thus, to obtain a high subregion matching quality and an accurate DIC calculation, 19 pixels will be the optimal value for R. Then, the size of the square subregion in this study is given by (2R + 1) × (2R + 1), namely, 39 × 39 pixels.

Parameter Setting of DIC Method
The DIC method, which is a full-field and noncontact-type deformation measurement method based on pattern matching, has been widely used in experimental mechanics [9,18,[32][33][34]. We also used this method to calculate the deformation field of the sample surface with the digital images of the original speckles of the sample surface obtained in the present experiment. The factors that influence the accuracy of the DIC measurement include the hardware system, the environment, the gray-level distribution in the image, the parameters set for calculation, etc. The systematic strain measurement error induced by self-heating of the cameras is remarkable and should be taken into account before measurement. According to a test of six cameras of four different types, a preheating period of at least 1-2 h is required to reduce the systematic strain measurement error [33]. In the present experiment, the preheating time was 3 h.

Influence of the Subregion Size on the Correlation Coefficient
A subregion centered at a point in the original image is chosen to calculate the correlation coefficient (CC) between the subregion and the moving subregion of the same size in the deformed image when the location of the point in the deformed image is searched. CC is an important parameter reflecting the subregion matching quality and the accuracy of the DIC calculation. CC ranges from 0 to 1 and increases as the subregion matching quality improves. It also depends on the subregion size ( Figure 3). The average CC attains a peak value when the radius of a subregion R is 19 pixels. Although the standard deviation of CC decreases as R increases, the decrease of the standard deviation of CC is remarkable when R is less than 19 pixels but is inconspicuous when R is greater than 19 pixels. On the other hand, because the computational time and space filtering effect increase as the size of a subregion increases, a smaller subregion results in a better DIC calculation under the premise of a desirable CC. Thus, to obtain a high subregion matching quality and an accurate DIC calculation, 19 pixels will be the optimal value for R. Then, the size of the square subregion in this study is given by (2R + 1) × (2R + 1), namely, 39 × 39 pixels. Uncertainties will be generated during subregion matching when the fault trace and subregions intersect. The distance between a measurement point and the fault, which makes the fault intersect the subregion centered at the measurement point, is called the fault influence length and is denoted by LA ( Figure 4a). As the fault is 45° oblique in the images, LA is half of the length of a subregion

Influence of the Subregion Size on the Fault Displacement Measurement
Uncertainties will be generated during subregion matching when the fault trace and subregions intersect. The distance between a measurement point and the fault, which makes the fault intersect the subregion centered at the measurement point, is called the fault influence length and is denoted by LA (Figure 4a). As the fault is 45 • oblique in the images, LA is half of the length of a subregion diagonal, namely, LA = √ 2(2R + 1)/2. In order to explore the influence of the intersection between the fault trace and the subregions on the fault displacement measurement, three subregion radii (9, 16, and 50 pixels) were used to calculate the fault displacement during a stick-slip instability via the DIC method. Figure 4b shows the change in the fault displacement versus the distance between the fault and the measurement point for the three values of R. In Figure 4b, the common feature of the three lines is that the fault displacement increases followed by a steady period as the distance between the fault and the measurement point increases. The distinct difference between the three lines is that the distance between the fault and measurement point at which the fault displacement begins to be steady is larger with increasing R. This distance is equal to LA. As the fault displacement should be equal at the same position during the same stick-slip instability, this difference cannot reflect the true displacement of the fault. This is caused by the intersection between the fault and the subregion within LA. Therefore, the distance between the fault and the measurement point should not be less than LA. Thus, the distance between the fault and the measurement point should become larger with increasing R. Doing so will eliminate the intersection between the fault trace and the subregions. On the other hand, to accurately obtain instances of slip along the fault surface, the selected points used for the fault displacement measurement should be located as close as possible to the fault. Accordingly, the subregion size must be small to ensure that these measurement points are as close as possible to the fault. Therefore, for a chosen subregion size of R, the optimum distance from the measurement point to the fault is LA. Thus, LA will be 27.6 pixels when R is 19 pixels. As the length of one pixel is equal to 149.5 µm, the length of LA will be 4.2 mm, which is the distance between the fault and the measurement lines denoted by l 1 and l 2 in Figure 1b. the fault trace and the subregions on the fault displacement measurement, three su bregion radii (9,16, and 50 pixels) were used to calculate the fault displacement during a stick-slip instability via the DIC method. Figure 4b shows the change in the fault displacement versus the distance between the fault and the measurement point for the three values of R. In Figure 4b, the common feature of the three lines is that the fault displacement increases followed by a steady period as the distance between the fault and the measurement point increases. The distinct difference between the three lines is that the distance between the fault and measurement point at which the fault displacement begins to be steady is larger with increasing R. This distance is equal to LA. As the fault displacement should be equal at the same position during the same stick-slip instability, this difference cannot reflect the true displacement of the fault. This is caused by the intersection between the fault and the subregion within LA. Therefore, the distance between the fault and the measurement point should not be less than LA. Thus, the distance between the fault and the measurement point should become larger with increasing R. Doing so will eliminate the intersection between the fault trace and the subregions. On the other hand, to accurately obtain instances of slip along the fault surface, the selected points used for the fault displacement measurement should be located as close as possible to the fault. Accordingly, the subregion size must be small to ensure that these measurement points are as close as possible to the fault. Therefore, for a chosen subregion size of R, the optimum distance from the measurement point to the fault is LA. Thus, LA will be 27.6 pixels when R is 19 pixels.
As the length of one pixel is equal to 149.5 μm, the length of LA will be 4.2 mm, which is the distance between the fault and the measurement lines denoted by l1 and l2 in Figure 1b.

Normalized Length of Preslip Zones
The locations of the fault where the dextral strike-slip component of the cumulative fault displacement is greater than 5 µm (the error of the fault displacement measurement) are considered as preslip zones, while other locations are considered as static zones. The length of a preslip or static zone can be calculated by the product of the number and spacing (0.2 mm) of the fault displacement sensors within the zone. By doing this, the total length of the preslip zones (L r ) can be obtained. The normalized length of preslip zones (L N ) [18] is defined as follows: where L is the length of the fault. The value of L N is between 0 and 1. The higher its value is, the greater the proportion of the length of the fault is occupied by the preslip zones. Thus, L N indicates the level of the expansion of the fault preslip zones.

Normalized Number of Isolated Preslip Zones
The total number of the fault displacement sensors is N (N = 1761), then the maximum number of the isolated preslip zones will be N m = N/2 (when N is an even number) or N m = (N + 1)/2 (when N is an odd number). N m means that the length of each preslip/static zone is the same and equal to the minimum value (0.2 mm) and the preslip and static zones are located alternatively. Therefore, the actual number of the isolated preslip zones (N r ) is not greater than N m . The normalized number of isolated preslip zones is defined as follows: From Equation (2), the value of N N is between 0 and 1, which represents the degree of dispersion of the preslip distribution along the fault. The smaller the N N is, the more concentrated the preslip zones distribute along the fault.

Normalized Information Entropy of Fault Displacement Direction
The normalized information entropy of fault displacement direction (S N ) is used to describe the dispersion degree of the displacement direction of each part of the fault based on the concept of Shannon entropy [35]. S N is defined as follows [18]: where Q is the total number of divided direction intervals within 360 • (direction interval of 5 • was used in this study and Q = 72 accordingly). The use of Q as the base of the logarithm makes S N normalized. p(i) is the probability that the direction of the fault displacement measured by the fault displacement sensors within the ith direction interval. p(i) is defined as follows [18]: where N i is the number of fault displacement sensors with fault displacement direction within the i-th direction interval; N is the total number of fault displacement sensors (N = 1761). Condition of p(i) = 0 will be omitted when Equation 3 is calculated. The value of S N is between 0 and 1. The smaller the S N is, the more concentrated the distribution of the fault displacement directions is.  Figures 5 and 6 show the spatiotemporal evolution of the preslip zones and changes of the normalized statistical parameters with time during the Events E1 and E2, respectively. The zero time represents the onset of the stick-slip instabilities. The yellow areas in Figures 5a and 6a represent the preslip zones, while the blue areas represent the static zones. The preslip zones expand independently prior to the onset of Event E1, but they migrate along the fault (the red dashed arrow represents the migrating direction) as they expand before the onset of Event E2. Figures 5b and 6b are the changes of the average cumulative fault slip (d), L N , N N , and S N with time during the Events E1 and E2, respectively. N N is displayed after magnified 10 times to ensure the change of it is visible. As shown in Figures 5b and 6b, before the onset of the stick-slip instability, d accelerates continuously, L N shows an increasing trend, N N shows a transitional trend from increase to decrease (the peak values of N N are −0.442 s and −0.864 s prior to Events E1 and E2, respectively), and S N shows a decreasing trend. The evolutions of L N and N N indicate that the expansion of the preslip zones along the fault begins with an increase in the number of isolated preslip zones but changes to the linkage of the isolated preslip zones. The peak value of N N indicates that linkage among the preslip zones begins to be a dominant factor controlling the evolution of preslip.    Thus, the initiation of linkage of preslip zones can be identified by the peak value of N N . The decrease in S N indicates that the displacement direction of each part of the fault gradually tends to be consistent. In addition, the changes of the normalized statistical parameters have the following common features: the rate of the increase of L N and the rates of the decline of N N and S N become slower as the instability approaches.

Fluctuation of the Normalized Statistical Parameters
The preslip zones exhibit a migrating propagation in Event E2 compared with Event E1. Accordingly, the evolutions of the normalized statistical parameters show more conspicuous fluctuations in Event E2 than that in Event E1 before the N N reaches its peak values. The fluctuations of the normalized statistical parameters appear as a pause or a reverse in their increase or decrease trends (Figure 6b). In addition, the three normalized statistical parameters exhibit different fluctuations. As shown in Figure 6b, L N and N N show a steady rise that are consistent with their general trends during the AB period, while S N appears to a pause that is inconsistent with its general decreasing trend. However, during the BC period, L N and N N show a pause or decline that are inconsistent with their general increasing trends, while S N shows a steady reduction that is consistent with its overall trend. These differences are contributed to the continuous acceleration of d and the double types of propagation of the preslip zones in Event E2. In the AB phase, the propagation of the preslip zones is dominated by the independent expansion. As a result, L N and N N increase due to the generation of new preslip zones. However, the directions of the fault displacement within the new preslip zones are dispersed due to the comparable values of the strike-slip and normal-slip components, which subsequently cause the pause in S N . In the BC phase, the propagation of preslip zones is dominated by migration, which does not enlarge the preslip range. Accordingly, L N and N N are unchanged or reduced slightly. However, the continuous accelerating d within the unchanged preslip range means the acceleration of preslip, which increases the dextral strike-slip component of the fault displacement in the preslip zones and accordingly reduces S N .

Migration of Localized Strain Bands Controls the Migration of Preslip Zones
The normal and shear strains along the fault were calculated with the measurement points in l 2 . The duration before the stick-slip instability was drawn in Figures 7 and 8 to ensure the visible of the process which will be too small to be invisible compared with the large change during the instability. The spatiotemporal evolutions of the normal and shear strains prior to Events E1 and E2 are shown in Figures 7 and 8, respectively. The localized strain bands can be identified in the two events. The difference is that the localized strain bands show parallel to the time axis in Event E1, which is consistent with the independent expansion of the preslip zones. However, they are oblique to the time axis in Event E2. The black dashed arrows in Figure 8 have the same spatiotemporal locations as that shown in Figure 6a, which indicates the spatiotemporal correlation between the migrations of the preslip zones and localized strain bands. Therefore, the migration of the preslip zones is caused by the migration of the localized strain bands that represent the dynamic adjustment of stress along the fault.

Discussion
Synergistic activity along faults was observed in field conditions [21]. In our previous study, the spatiotemporal evolution of the fault displacement in the middle part of the fault was observed and the decrease of S N was revealed, which indicated the synergistic activity of the fault displacement before the instability [18]. The similar process was again observed in this study with fault displacement observations covering the entire fault, which showed the evidence for synergistic activity along the entire fault prior to its instability. The use of the evolutions of L N and N N revealed that the linkage of the isolated preslip zones is the mechanism of the transition from multiple preslip zones to single preslip zone prior to the stick-slip instability. The linkage of the preslip zones initiates at the point where N N reaches its peak value and enhances the connection of different sliding segments along the fault. As a result, the activities of different parts of the fault tends to be consistent and the spatial complexity of preslip reduces, which may be the mechanism of the synergistic activities of faults observed in field or laboratory conditions. The density of spatial distribution of sensors plays crucial role in uncovering the evolution of preslip in the laboratory. The use of one displacement gage mounted near a simulated fault could detect stick-slip preceded by preslip, which could not be readily observed at the specimen endcaps [2]. However, the use of multiple strain/slip gages can detect the propagation of a single preslip zone [3][4][5]. In recent years, multiple preslip zones were detected by using multiple types of gages (piezoelectric sensors, strain gages, and slip sensors) mounted on both the top and bottom sample surfaces [8], dense distribution of strain gages [20], or very dense distribution of fault displacement sensors via the DIC method [9,18]. The very dense distribution of the fault displacement sensors via the DIC method in this study again enable to obtain a detailed spatial distribution of the preslip zones that show a transition from randomness to concentration due to the linkage of them. These indicate the dependence of how much we know about the evolution of preslip on the development of the observation technology or methods.
Experimental studies have observed the interaction between the preslip and asperities, which revealed that the asperities can obstruct and promote the preslip before and after their failure, respectively [8,36]. The stability of the fault is actually supported by a few strong asperities [37,38]. The successive failure of these asperities causes a large amount of stress to accumulate on the last static asperity, which means that more strain energy needs to be accumulated to break the asperity. Therefore, the last static asperity will have an obstructive effect on the preslip behavior before its failure, which leads to the slower rates of increase of L N and decreases of N N and S N prior to the stick-slip instability as shown in Figures 5b and 6b.
The migration of localized strain bands along the fault (as shown in Figure 8) was similar in other experiments previously reported in our conference paper [39], the origin of which is still not clear (it may be caused by a complicated interactions between preslip and deformation of asperities, or the additional load generated by the press system under certain conditions). However, its effect on the preslip is obvious. The spatiotemporal correlation between the migration of the localized strain bands and the preslip zones shown in Figures 6a and 8 indicates that the dynamic adjustment of the stress on the fault controls the spatiotemporal variation of the preslip zones and leads to the fluctuations of L N , N N , and S N shown in Figure 6b. The fluctuations will probably bring uncertainty to the real-time identification of the peak N N which indicates the approach of the impending instability. Therefore, future study on the origin of the dynamic stress adjustment will be of great significance for reducing the uncertainty. Furthermore, such a future study may also implement the new method of natural time analysis which has been applied to laboratory measurements [40] and dynamical models [41,42] associated with stick-slip processes as well as to real seismic data of Japan [43,44], California [45], etc., in order to identify the approach of the system to the critical point. Notably, natural time has been recently used as a basis for a new method of estimation of seismic risk by Turcotte and coworkers [46][47][48][49].

Conclusions
The detailed spatiotemporal evolution of preslip zones was observed in this study, based on which three normalized statistical parameters were defined to quantitatively describe the evolution of preslip. The increase in L N and the transition from increase to decrease in N N indicate that the expansion of the preslip zones along the fault begins with an increase in the number of isolated preslip zones but changes to the linkage of the isolated preslip zones. The peak N N is an identifier of the initiation of linkage among preslip zones. The linkage of the preslip zones and the decrease in S N indicate the reduction of the spatial complexity of preslip, which is subjected to fluctuation by the migration of the localized strain bands controlling the migration of the preslip zones and is slowed down by the interactions between the preslip and asperities as the instability approaches.