Frictional anisotropy of 3D-printed fault surfaces

The surface morphology of faults controls the spatial anisotropy of their frictional properties, and hence their mechanical stability. Such anisotropy is only rarely studied in seismology models of fault slip, although it might be paramount to understand the seismic rupture in particular areas, notably where slip occurs in a direction different from that of the main striations of the fault. To quantify how the anisotropy of fault surfaces affects the friction coefficient during sliding, we sheared synthetic fault planes made of plaster of Paris. These fault planes were produced by 3D-printing real striated fault surfaces whose 3D roughness was measured in the field at spatial scales from millimeters to meters. Here, we show how the 3D-printing technology can help for the study of frictional slip. Results show that fault anisotropy controls the coefficient of static friction, with the friction coefficient along the striations being three to four times smaller than the friction coefficient along the direction perpendicular to the striations. This is true both at the meter and the millimeter scales. The anisotropy in friction and the average coefficient of static friction are also shown to decrease with the normal stress applied to the faults, as a result of the increased surface wear under increased loading.


I. INTRODUCTION
Active faults in the Earth's crust are complex systems along which earthquakes nucleate and propagate [1]. Faults hold structures and heterogeneities at all scales [2][3][4]. While they are often simplified to their simplest two-dimensional description (i.e., the fault plane), increasing complexity is now added to faults models [5]. It is indeed considered that, to fully understand seismicity in various areas [6][7][8][9][10][11][12][13], it is paramount to account for some disorder in the faults frictional properties such as secondary faulting, off-fault damage or roughness of the fault plane. Another degree of complexity is sometimes considered when modelling geological contacts and fault slip: the fault roughness and the possible anisotropy in their frictional properties. Morphological anisotropy is a known feature of faults, notably impacting the seismic waves velocity in their vicinity [14][15][16] or the mobility of natural and injected fluids [17] in the subsurface. Frictional anisotropy, interestingly, is also regularly studied in other fields than seismology, for instance the tribology of rubber tires [18,19], the strength of advanced adhesives [20], or the mitigation of water condensation [21]. It is also considered to play a major role in nature [22], for instance in the motion of numerous animals [20,23,24] and the hydration of some plants [25,26]. In most cases, frictional anisotropy derives from the existence of preferential topographical orientations on, at least, one of the contact surfaces [27,28]. The length scale for such structural directivity can be as small as micrometer [29] to nanometer [30,31]. In seismic faults, such preferential directions in their topography are observed at all scales [2][3][4] and originate from several processes. At the molecular level, rock forming crystals may display some frictional anisotropy. It is notably the case for antigorite [32], a mineral abundant in the Earth's upper mantle. At the mesoscopic scale, the shear strength of foliated rocks is known to be anisotropic [33,34], due to the oriented planes in their constitutive mineralogy. Fault zones in sedimentary basins are initiated by early fractures that often propagate in highly layered sediments. It can result [35] in an anisotropic ramp-flat morphology of these fracture surfaces. For more mature faults having accumulated enough displacement, and above a given length scale [36], the topography of the fault planes is also marked by slip induced wear, with striations and grooves of various wavelengths and amplitudes [37,38] oriented along the main direction of slip. If such morphological anisotropy of fault surfaces is well-known, its effect on the anisotropy of the frictional properties remains to be characterised. Such a characterisation of frictional anisotropy could also be of interest for other types of rock contacts than strictly seismic faults, in particular for shallow rock joints and fractures [39][40][41][42], whose three-dimensional geometry is key in geotechnical engineering and for the structural stability of many man-made constructions. Here, we study how the morphology of faults controls the static coefficient of friction and the anisotropy of friction with regards to the main stress direction during slip. To reach this goal, we produce 3D-prints of actual faults surfaces whose topography was measured in the field [43]. We perform friction experiments with plaster of Paris casts of these 3D-printed faults. Results show that the coefficient of static friction along faults is highly anisotropic, a property that should henceforward be considered in numerical models of slip on seismic faults. We also show that this anisotropy is stress dependent, and should decrease with depth (e.g., [44]).

II. 3D PRINTING AND PLASTER CASTING OF FAULT PLANES
The actual morphology of natural faults can be difficult to assess, even if their long wavelength structures can be inferred by surface or subsurface imaging techniques [38,45,46]. Yet, some fault planes are accessible to direct, high resolution, measurements, notably as they were exhumed by erosion and tectonic processes. For this study, we have used a series of digital fault surfaces. These fault roughness data were acquired with Light Detection and Ranging (LIDAR), laboratory laser profilometry, or white light interferometry [3]. These data are available on an online public database [43]. Should the reader hold some similar data, these authors welcome additions to this database. We have specifically selected fault roughness measurements performed on the Corona Heights fault [11] that outcrops near the Peixotto playground in San Francisco, California. These data cover surface areas with spatial scales in the range of millimeters to meters. Figures 1 and 2 show the fault surface at two spatial scales, one surface at the meter scale, defined on a 5 mm×5 mm grid, and one surface at the millimeter scale, defined on a 2 µm×2 µm grid. We will further on refer to these two surfaces as respectively S m and S mm . Already, one can notice some preferential directions in these topographies, and that the amplitude of fault roughness is, relatively to their size, somewhat larger at smaller scales (S mm ) than at larger scales (S m ) [47]. For our tests, we chose to limit these anisotropic surfaces to a circular sample geometry. We also applied a mild running-window median filter to smooth out spikes in the measured surfaces that could be associated to measurement noise. The window length of the filter was 10 space steps, accounting for 5 cm for S m and for 20 µm for S mm . In order to run the friction experiments, we generated some opposing surfaces to the ones presented in Figs. 1 and 2. These opposing surfaces could not be measured, as the actual fault walls that were facing S m and S mm are now eroded. To reconstruct them, we have applied the following transformation to the 3D coordinates (X, Y, Z) of S m and S mm : where X , Y and Z are the coordinates of the generated opposing surfaces and (X, Y ) give the map location of a given surface point of elevation Z, as represented  [43]. This surface, called Sm, has a radius of 1 m and is defined on a 5 mm grid with a 1.25 mm elevation resolution. A parametric angle θ is defined from the main groove direction.
FIG. 2. White light interferometry measurement of the topography of the Corona Heights fault plane at the millimeter scale. This surface, called Smm, has a radius of 1.5 mm and is defined on a 2 µm grid with a 0.025 µm elevation resolution. A parametric angle θ is defined from the main groove direction.
in Fig. 1. We have thus assumed that the missing fault walls are complementary to the measured ones, so that, when pressed together before the friction tests, they form a bulk with negligible aperture between the two blocks. Such assumption for natural faults would only be partly verified. When having accumulated enough slip, a gran- ular layer of gouge material may there have formed, and the two opposing sliding surfaces may not always perfectly match. However, our assumption is relevant for the youngest faults with a small amount of slip. We have also assumed that erosion did not significantly alter the fault plane, such that the measured topography is representative of the one of an actual buried fault. For the Corona Heights fault, this assumption is valid because the fault offsets silica-rich chert rocks with a high resistance to weathering. After having obtained the surfaces, we isotropically (i.e., with the same factor in all directions) down-or up-scaled S m and S mm to fit a standard 4 cm diameter disk that matches the clamp size of our shear deformation apparatus. We also re-gridded the surfaces to match the lateral resolution of our 3D-printer ('Ultimaker 2 Extended+' [49]) that has a nozzle size of 250 µm. The four surfaces (two fault surfaces and two opposing surfaces) were then 3D-printed into polylactic acid (PLA) material, as shown in Fig. 3. It should be noted that, even when designed to be flat, printed objects can present a natural roughness [50], at a scale however smaller that the grooves observed on the printed faults. These intrinsic imperfections shall be comparable to 60 µm, the elementary thickness of the PLA layers deposited by our 3D printer. In comparison, the 2D standard deviation of the elevation in our printed objects topography are 0.66 mm for S m and 1.7 mm for S mm . The maximal elevation of these objects are, respectively, 3.7 mm and 8.3 mm. We thus consider that the small scale roughness (∼ 60 µm) from the printer's limit in resolution has a second order effect on the frictional properties of the surfaces. Although we could have performed the friction experiments with the plastic pieces produced with the 3D printer, we have rather produced samples of plaster of Paris (gypsum) blocks moulded from the plastic faults. Plaster is known [51] to be a reasonable model of porous brittle materials, and the main goal of this casts was to work with a rock-like material, notably because it may wear and deform differently than plastic under shear. Their fragile nature, and the potential friction-induced wear that the plaster was subjected to in our experiments, made us use new casts for each experimental realisation. The casts were generated with the following protocol: five volumes of water and eight volumes of powder of plaster of Paris were mixed and poured over the plastic moulds, then let to dry during one and half hour. The moulds, an example of which is shown in Fig. 3, were sprayed before each cast with a thin layer of silicon grease  (3) necessary to obtain a null normal loading when required (not mounted on the picture); and two horizontal and vertical sliders (4) used to force the motion in the direction of interest while allowing for vertical displacement. The shear force FT is applied on the top fault wall in the e θ direction, while the bottom surface is kept fixed. It is measured by a Sauter ® FH500N force gauge [48]. The normal force FN is applied by a dead weight that acts oppositely to the e Z direction) on the top surface. For a visualisation purpose, the two plaster blocks (1 and 2) are offset vertically in the schematic, whereas in the experiments these two surfaces matched, with a quasi-null aperture.
to avoid some of the fine plaster details to stick to the plastic during the mould release. The last step in the casts preparation was to dry them in an oven at a temperature of 40 • C for one hour. As a result, we produced fault planes in blocks made of plaster of Paris, as shown in Fig. 4.

III. EXPERIMENTAL SET-UP AND EXPERIMENTAL CONDITIONS
The shear apparatus used to perform the friction tests is shown in Fig. 5. The two complementary surfaces are pressed together and mounted one on top of the other between the clamps of the shear apparatus. A normal force F N is applied on the top surface by using adjustable weights. In addition, a spring system of stiffness 625 N m -1 allows, if desired, to compensate for the machine empty weight of 13.7 N (i.e., the normal weight transmitted to the friction surfaces by the machine top clamp and structure and the top cast when no extra mass is used). A tangential driving shear force F T is then applied to the top fault wall in a given direction of the (X, Y ) plane. The amplitude of the force is measured by a Sauter ® [48] force gouge. The shear direction direction is defined by the angle θ ∈ [0 2π] from the direction of main grooves on the fault surface, as defined in Figs. 1 and 2. An horizontal mechanical slider makes sure that the friction is evaluated in the direction of interest only, and a vertical slider allows upward or downward displacement of the top surface. While these sliders would ideally be perfectly lubricated, we have estimated their frictional resistance at the sliding velocity of our experiments, F c = 4 N±0.5 N by performing a friction test with no fault installed in the machine (that is, with only air between the two clamps represented in Fig. 5). We characterise the anisotropy of the static friction coefficient for both the S m and S mm surfaces by performing a series of experiments where we vary the direction of loading with respect to the grooves. We ran friction tests every 30 • on both plaster faults. At each angle, the experiment was repeated at least three times (with new casts) in order to ensure that results are reproducible. The standard deviation computed on these multiple measurements (typically 5 to 10 N) was used to compute the error bars on our characterised coefficients of static friction. S m was sheared under a normal stress σ N = 10.9 kPa, while the tests performed on the rougher surface S mm were performed under σ N = 26.5 kPa. A total of 76 experiments were performed, 37 using S m and 39 using S mm . At the onset of slip, the static friction coefficient is defined using a standard Coulomb's law [52]: where µ s is the coefficient of static friction and max [F T ] is the maximum tangential force applied at the onset of slip. In the following, we will also consider the mean driving and normal stresses, denoted σ T = F T /(πr 2 ) and σ N = F N /(πr 2 ), where r = 2 cm is the radius of the cast. Figure 6 shows a typical measurement of a friction test, from which max [F T ], and hence µ s , are calculated. For the record, the target speed of the test bench speed (that is, the demanded slip velocity) was fixed to a constant 1.3 mm s -1 .

A. Friction anisotropy
The results are presented in Fig. 7. On both fault samples, one can observe the strong anisotropy of the coefficient of static friction, with the maximum value of σ S being about four times higher than its minimum for S m and about three times bigger for S mm . In most experiments, the minimal friction is obtained along the main groove direction (i.e., at θ = 0 • or θ = 180 • ), and the resistance to shear is higher perpendicularly to this direction. The maximum of value is however never obtained exactly at θ = 90 • or θ = 270 • but rather along a neighbouring direction. Local maxima are indeed obtained for θ = 60 • or θ = 240 • when shearing S m and for θ = 120 • or θ = 300 • when shearing S mm .

B. Damage and stress dependence
Most of our experiments were destructive, with visible wear on the plaster samples after the shearing tests. The observed damage consists either in the formation of plaster powder (gouge) or in the rupture of topographic highs of the faults surfaces. Part of it might have initiated at the onset of the fault displacement (and hence be related to the static friction), while some of it has rather been induced by the subsequent dynamical friction. Figure 8 shows some examples of these damage types. Wear of seismic faults has been studied (e.g., [53,54]), because this process may lubricate faults during slip (e.g., [55,56]). The present study focuses on the measurement of the coefficient of static friction and on its anisotropy, but we suggest that our 3D-print-based setup could also enable the quantitative characterisation of damage during sliding along analogue fault surfaces. We here keep to a qualitative assessment of which parts of the surfaces were mainly worn during each experiment. It seems that most of the shear resistance of the Corona Heights fault, at the millimeter scale (S mm ), arises from its grooves. By contrast, the friction of S m (representing a metric scale) is dominated by one chip in its morphology (i.e., the main topographic low in Fig. 1). This chip being on the edge of the 3D-printed surface, a finite size effect is certainly at play in the results reported in Fig. 7 (top), notably explaining the strong asymmetry in µ S for the opposite directions θ = 60 • and θ = 240 • . Because the overall friction is likely to be affected by the surface wear, and because this wear is likely stress dependent, we have performed some friction tests on S m under various loads σ N , at a given angle θ = 90 • . The results are shown in Fig. 9. At the highest tested stresses, µ S seems, to an extent stress independent, with its mean variations lying within the measured error bars. While this result is compatible with the classical Coulomb theory (e.g., [52]), one can observe, over a wider range of normal stresses, a consistent decrease in the friction coefficient with a higher normal stress. It could, in part, emanate from some limitations in our experimental setup. For instance, at lower σ N , the internal friction of the device (F c ) accounts for a more significant portion of the total measured tangential force, and the friction characterisation could thus be less accurate. We have however run similar tests on flat plaster surfaces, showing no significant variations in µ S over the same stress range (see Fig. 9). The drop in friction coefficient with σ N is then likely related to the increased damage under higher normal loading (see Fig. 8), reducing the overall shear resistance as asperities are easier to break. Such a drop of µ S with the normal stress has already been reported for rough contacts [57], as the static coefficient of friction is not an absolute material constant [57,58]. We have then assessed the effect of the normal stress σ N on the anisotropy of the static friction coefficient. We performed frictional tests on the four poles of S m (θ = 0 • , 90 • , 180 • and 270 • ) at a higher load (σ N = 171 kPa) than what was used before (i.e., σ N = 10.9 kPa, as reported in Fig. 7). The newly measured coefficients of static friction are shown in Fig. 10. One can notice the reduced friction anisotropy at high σ N . The ratio between the maximum and the the minimum value of µ S indeed drops from 3.1 at σ N = 10.9 kPa to 1.5 at σ N = 171 kPa. This result suggests that the frictional anisotropy of faults is smaller at depth. This concept is naturally linked with the seminal Byerlee [44]'s law, stating that the roughness of fault planes (as well as the type of their constitutive rocks) has less effect on their frictional properties at larger depths.

V. DISCUSSION AND CONCLUSION
Here, we show how the multiscale anisotropy of fault plane topography leads to an anisotropy in the frictional properties. Results confirm that seismic faults are prone to slide along some preferential directions. The direction that is the most likely is the one that faults have previously slid along, and which has shaped some guiding grooves in their morphology. Yet, displacements following other orientations are possible. Predicting the rupture direction of the next earthquake on a fault is thus not only dependent on assessing the main regional stress. The question should rather be along which orientation a rupture criterion [52] will first be exceeded. Such a subtlety might be of little importance for mature faults for which the stress principle directions have not changed with time, because, in this case, the main stress is likely to act along the lowest coefficient of friction anyway. Yet, it could be paramount for faults under a changing geological load, where this alignment is not verified, or for immature faults, where the slip could be mainly governed by the anisotropy of early surfaces (i.e., where the slip does not coincide with the stress principal direction, but is non-associated). Earthquakes occurring along abnormal directions (i.e., not in agreement with the local stress state) do occur [59][60][61], and their understanding might be eased by accounting for the possible frictional anisotropy of their surfaces [62]. While we have here only measured the static coefficient of friction, we note that the same considerations could be applied to characterise the coefficient of dynamic friction. Hence, not only the initial slip direction of an earthquake could be impacted by frictional anisotropy, but its complete trajectory [28]. We have, additionally, measured how the anisotropy in friction becomes less significant when the normal stress acting on a fault increases (i.e., with the fault depth), in general agreement with Byerlee's law [44]. Such an effect likely derives from the stress related changes in rupture rheology and in damage type. The transition from a highly anisotropic to a relatively isotropic regime should typically occur when local stresses on the fault reach the yield stress of the material, σ y ∼ 5 MPa [51] in the case of plaster. This is about two orders of magnitude above the transition σ N at which we observed a strong reduction in anisotropy (see Fig. 10), but our simple computation of σ N does not account for the potential strong stress concentrations at play in our faults. Considering that the strength of rocks (e.g., [63]) is about two orders of magnitude larger than that of plaster, fault frictional anisotropy could thus be only at play at pressures less than about σ N ∼ 10 MPa. This would correspond to the shallowest faults, at depths less than σ N /(gρ) ∼ 500 m, where ρ is the volumetric mass of rocks (∼ 2000 kg m -3 ) and g the gravity acceleration. Some care should however be taken when deriving such a conclusion by analysing resized samples as ours (40 mm diameter samples representing meter or millimeter topographies), as the way matter breaks is length-scale dependent (e.g., [47]). Additionally to the assessment of the stability of shallow seismic faults, the characterisation of the frictional anisotropy of rock surfaces may be of importance in geotechnical engineering, for instance, for the stability of tunnels and foundations. There, the intrinsic strength anisotropy of foliated rocks is well studied [33,34]. Our work shows how one can also characterise the mechanical anisotropy of rough rock contacts, for instance, along joints [39][40][41][42] and fractures [64,65] between or inside rock formations. A main point of this manuscript was, finally, to illustrate how the 3D-printing technology can help with new experimental designs in Earth Sciences, and this technology is getting a growing attention from the community [66][67][68][69]. A direct continuation of the present work, for instance, could be to 3D-print and to test some faults surfaces beforehand filtered with various band-pass filters, in order to understand how the various wavelengths of the topography contribute to the global static friction coefficient, to the dynamical friction coefficient and to analyse the spatial distribution of the fault wear produced under various stresses and amounts of slip.