Scalar model for frictional precursors dynamics

Recent experiments indicate that frictional sliding occurs by nucleation of detachment fronts at the contact interface that may appear well before the onset of global sliding. This intriguing precursory activity is not accounted for by traditional friction theories but is extremely important for friction dominated geophysical phenomena as earthquakes, landslides or avalanches. Here we simulate the onset of slip of a three dimensional elastic body resting on a surface and show that experimentally observed frictional precursors depend in a complex non-universal way on the sample geometry and loading conditions. Our model satisfies Archard's law and Amontons' first and second laws, reproducing with remarkable precision the real contact area dynamics, the precursors' envelope dynamics prior to sliding, and the normal and shear internal stress distributions close to the interfacial surface. Moreover, it allows to assess which features can be attributed to the elastic equilibrium, and which are attributed to the out-of-equilibrium dynamics, suggesting that precursory activity is an intrinsically quasi-static physical process. A direct calculation of the evolution of the Coulomb stress before and during precursors nucleation shows large variations across the sample, explaining why earthquake forecasting methods based only on accumulated slip and Coulomb stress monitoring are often ineffective.


I. Introduction
The classical laws of friction, due to Amontons and Coulomb, postulate that a body resting on a surface can be displaced only by applying a shear force larger than a static friction force, which is proportional to the normal load and independent of the apparent area of contact. Recent research has challenged this understanding of friction, showing that macroscopic slip is due to the formation and propagation of detachment fronts through the contact interface [1]. The nature of these fronts and their speed depend on the way shear is applied to the sample and on its geometry [2,3], a particularly compelling issue in view of the long held assumption of independence of friction on the sample shape and size. It is particularly intriguing that, in some cases, localized sliding precursors nucleate long before the applied force reaches the static friction force at which the front propagates through the entire contact interface [1, 4,5]. Numerical simulations of friction models in one [6][7][8][9][10] and two dimensions [11][12][13] allow to study the main features of the spatio-temporal dynamics of precursors. These numerical works have mainly focused on the qualitative dynamical aspects of propagation, reproducing the different dynamical regimes observed in experiments [6][7][8][9] and the nucleation of the fronts under various loading conditions [11].
Based on the results of experiments [4] and numerical simulations [11] it was suggested that frictional precursors evolve according to universal laws: the sample size and normal load dependences of precursors lengths can be rescaled away and different experiments can be collapsed a single master curve. Establishing universal forms for slip precursors would be particularly important for earthquake forecasting [14]. Slip or stress accumulation on faults has been often observed to accelerate close to large earthquakes [15][16][17], but detailed predictions based on this are considered to be unreliable [18,19]. It is therefore extremely important to better clarify the conditions leading to precursors and confirm their universality.
Another puzzling aspect revealed by experiments is an apparent violation of the Amontons-Coulomb laws: direct measurement of shear τ and normal stresses σ close to the frictional interface indicated regions where the Coulomb stress τ C = |τ | − µ |σ| is positive without inducing detachment [2,3]. This result suggests that the friction coefficient µ might not be a well defined material constant as conventionally assumed. Scalar models are commonly used to study the planar crack front propagation in disodered elastic media [20,21], in quasi-two dimensional geometries [22] and under antiplane shearing conditions [23]. On the other hand, recent experiments have been provided the evidence that classical shear cracks singular solutions, originally devised to account for brittle fractures, offer a quantitative excellent description of the static-to-dynamic friction transition [24].
Here, combining the solution of three dimensional scalar elastic equations in finite geometries with simple contact mechanical rules for local slip at the frictional interface, we reproduce accurately the experimentally observed evolution of the contact area as the sample is loaded.
In this way, we obtain a complete picture of the role played by sample geometry and loading conditions on the precursors nucleation. Moreover we show that precursors originate from stress gradients on the contact interface and are therefore absent when loading is applied uniformly through the top of the slider. Disorder induced precursors nucleation of the kind predicted for sliding thin films and monolayers [25] should be strongly suppressed in three dimensions due to long-range elastic interactions which make the coherence length extremely large [26]. When shear is applied on the slider side, however, we observe precursors whose evolution depends in a non-universal way on the sample geometry. The occurrence of universal profiles is explained by the symmetries of the interfacial shear stress obtained analytically.
Despite its quasi-static nature, our model incorporates the Achard's law and the Amontons' first and second laws, reproducing several key features observed in experiments including the discrete stress drops observed in correspondence of the slip precursors. Most importantly, our model can help to assess which experimental feature can be attributed to the static elastic equilibrium, and which instead is a pure dynamical out-of-equilibrium aspect.
Our calculations reproduce the experimental interfacial stress profiles detected at the frictional interface, before slip. We discuss the large fluctuations of the internal stresses during precursors activity in the bulk of the material, and we provide the numerical and analytical evidence of this large heterogeinity. This observation, substantiated also by finite element model simulations, suggests that drawing firm conclusions based on the value of the Coulomb stress measured away from the contact interface, both in laboratory experiments and in earthquake faults, could be problematic.

A. Simulations for different sample sizes and loading conditions
Following Ref. [4], we first study how precursors depend on L x and on the normal load F N when F lat S is applied through a rod placed on the trailing edge, at height h = 6mm.
Experimental evidence suggests that the precursor size obtained for different values of L x and F N can be collapsed into a single master curve when normalized by L x and plotted versus F lat S /µF N . Our numerical results reproduce quantitatively the experimental findings as shown in Fig. 5(a). In our model, however, we are able to change L x over a wider range than in the experiments, revealing that data collapse is in fact only approximate (see inset of constant the other parameters: in all cases front precursors exhibit a dependence on the sample dimensions. We have also changed L x and L y holding their ratio L x /L y unchanged ( Fig. S6), L x and L z with L x /L z constant (Fig. S7), or L y and L z with L y /L z constant ( Fig. S8). Again, data collapse is not obtained, indicating that for this loading condition the precursor lengt depends in a non-trivial way on the sample dimensions (L x , L y , L z ). The general trend however is that precursory activity tends to decrease as the varying dimension is increased: for a larger sample we typically need a larger shear force to observe a precursor of a given length.
Experimental results in Ref. [4] also suggest that the height h at which the lateral force is applied to the sample trailing edge has no influence on the evolution of the front precursors. While this is true for the range of h used in experiments (see Fig 5(b)), when we increase h further the curves no longer collapse. In particular, we find that the lateral force needed to nucleate the first precursor increases considerably with h (see the inset of Fig   5(b)). Remarkably, this effect persists when we increase both h and L z leaving their ratio constant (Fig. S9). Yet, experiments have provided the evidence that the precursors length advances by periodical discrete leaps of roughly equal size, which take place at nearly constant increments of F S . Moreover, this periodicity exhibits an apparent scaling with h, becoming larger with increasing h [27]. While the envelope of the curves reproduced by our model do show periodicity in the increments of the precursors sizes, the size of these discrete jumps and the corresponding increments of F S seem to remain unaffected by varying h, at 5 least within the range of heights used in the experiments.
We explore further the dependence of precursors on the sample geometry by considering a different loading condition in which the lateral force is applied uniformly on the sample side surface (2∆h = L z ). In this case, we find that the precursors are size independent when we vary L x and L y keeping their ratio L x /L y constant ( Fig. 6(a)) or L x and L z with constant L x /L z (Fig. 6(b)). When we vary instead L y and L z , keeping constant L y /L z , no universality is found and precursors again tend to disappear for large sample sizes ( Fig.   6(c)). A similar effect is obtained using mixed mode loading as in Refs. [2,3], applying simultaneously a shear force on the top surface and on the trailing edge. As the ratio between both forces n ≡ F top S /F lat S increases, the length of the precursor shrinks ( Fig. 6(d)) and disappear when loading is only applied on the sample top plate.
Precursors are defined by detecting the decay of the real area of contact. This feature is perfectly reproduced by our model, and, roughly speaking, it is ultimately due to the detachment of regions of the frictional interface satisfying the local static friction condition (19). Thus, the dependence of the precursor envelope profile on the sample geometry should reflect the properties of the Coulomb stress across the entire contact interface. In the Supplementary Information (sec. I), we discuss how some general aspects of the precursor shape can be deduced from the symmetry of the Green function.

B. Normal and shear stresses at the interface and the Amontons law
Direct experimental measurements of shear and normal stress profiles close to the contact interface show that the Coulomb stress can exceed zero locally, without inducing any detachment front, precursor or local slip [2,3]. This is puzzling since it would represent a local violation of the Amontons-Coulomb law, suggesting that the friction coefficient might not be a material constant. In our model, however, the local and global friction coefficient µ is fixed across the whole interface, and local detachment occurs if τ C (x, y, 0) > 0 by construction (Eqs. (19) and (I5)). Yet, this apparent contradiction can be resolved by noting that local stresses in Refs. [2,3] are measured on a reference plane located at a height of z P = 2mm above the frictional interface.
Thanks to the analytical solvability of our model we can compute the shear and normal stresses at any points (x, y, z) of the slider bulk: this is perfomed in the Suppl. Mat. sec. H 6 (see Eqs.(H2), (H12) or Eqs.(H17),(H20)). Calculating the stresses on the plane z P = 2mm yields a good quantitative agreement with experiments (Fig 7(a), S10). In particular, the curves shown in Fig. 7 represent the shear and normal stresses averaged over the y direction (τ (x, z P ) = Ly 0 dy τ (x, y, z P )/L y and σ(x, z P ) = Ly 0 dy σ(x, y, z P )/L y ), along the entire sample 0 < x < L x , just before the onset of the first precursor, i.e. when no detachment is yet present at the contact interface. This can be seen also from Fig.S10 where the full quasi-static dynamics of τ (x, z P ) and σ(x, z P ) are plotted. The corresponding Coulomb stress on the same plane (τ C (x, Eq.(H22)) prior to the nucleation of the first precursor event is reported in Fig.7(b) (solid green line) showing again a good agreement with the experimental data.
Our result may suggest the observed apparent violation of the Amontons first law [2,3] could be due to the fluctuation undergone by the internal stresses in the material bulk, even in the vicinity (but not on) the slider frictional interface. Defining a local friction coefficient as the ratio µ(x, y, z) = |τ (x,y,z)| |σ(x,y,z)| , is not an eligible procedure if the point (x, y, z) does not lie on the frictional plane (x, y, 0). To substantiate this statement, in Fig.7(b) we compare the y-averaged Coulomb stress on the plane z P = 2mm above the slider-bottom surface (τ C (x, z P ), solid green curve) with the corresponding quantity at frictional the interface τ C (x) = Ly 0 dy τ C (x, y)/L y ( solid magenta curve). As it can be clearly seen, the Coulomb stress value may suffer large fluctuations according to the sample position where it is measured. Although the authors of the experiments in [2] were careful to perform the measurements at locations x to "avoid the effects of large stress gradients", the agreement shown in Fig.7 and the analytical calculations in Suppl. Mat. sec. H demonstrate that large fluctuations in the (positive) Coulomb stress appear in the bulk of the material, and are mainly due to the internal shear stress gradient resulting from the lateral shear applied. A pictorial intuitive illustration of this argument is reported in Fig. 8, where the y-averaged Coulomb stress τ C (x, z) is shown for several values of the adiabatic force F S , see also Movie S2. Regions where τ C (x, z) exceeds zero occur even before precursors nucleation ( Fig. 8(a)) and even far from the trailing edge. This may lead one to believe incorrectly that the friction coefficient is not a material constant. In the same Fig.8 (bottom panels) we report the full Coulomb stress τ C (x, y) at the plane (z = z P ), showing that the average over y is indeed a correct protocol to average out the noise-induced fluctuations.
The bulk fluctuations of Coulomb stress are also present while considering a fully tensorial 7 elasticity model. As a matter of fact, in the Suppl. Mat. sec. J we report shear and normal stress calculated by means of a finite element model (FEM). In Fig.S12 we plot τ C (x, y = L y /2, z) arising from FEM simulations: large gradients of Coulomb stress make any claim about the local value of µ highly questionable. This is seen before any detachment occurred at the interface (panel (a)), and also when a portion of the contact area is disconnected from the rough surface (panels (b)-(d)). We notice that in this case no analytical calculations can be carried out in the fashion of Suppl. Mat. sec. H, but the shear and normal stresses are numerically obtained by means of the FEM software. In particular we stress once again that the the normal component of the positive Coulomb stress σ(x, y, z) is largely influenced by the shear force F S as opposite to Eq.H12.

C. Role of disorder in front nucleation
The role of the substrate roughness heterogeneity on the nucleation of front appears, from Eqs.(I1),(I3), rather complex. In general one can say that heterogeneity amplifies and modulate the different contributions to normal and shear interface stresses σ surf and τ surf . For narrow roughness distributions (in most of the experiments the roughness appears to be a very well-controlled parameter), we do not expect that the substrate disorder may play a major role promoting or suppressing the precursors dynamics, at least not comparable to the role expressed by the force-induced stress gradients. This is what clearly appears from our analysis, in step with the simplified 1D model [10] and with the experimental outcomes. It is possible, however, that large roughness fluctuations may induce internal stress gradients leading to precursors nucleation, even in regions far from the trailing edge. This is a particularly interesting issue, since the precursor could originate as a stable detachment droplet, irrespective of the type of loading exerted (whether lateral or top shearing). The key question is therefore to determine the conditions for which a detachment droplet constitutes a meta-stable state. The question could be addressed by defining the contact interface energy density as [28], and the energy change associated to the transition from an initial stable configuration to a second on which the droplet has formed: ∆E(r) = Σ(r) dxdy τ surf (x, y)u x (x, y) + σ surf (x, y) u , where Σ(r) represents the contact surface configuration including the nucleated droplet of average size r. A detachment 8 droplet configuration will be stable if the energy penalty ∆E(r) has a positive maximum for some r. Now, asking which physical conditions allow the formation of a stable droplet, means which sample dimensions L x , L y , L z , roughness w 2 and force F S (< µF N ) give a ∆E(r) with a positive maximum (with r < L x , L y , L z ). Unfortunately, due to the intricacy of expressions (I1) and (I3), we could obtain the answer only by numerical simulations.
However, albeit one cannot completely exclude that detachment regions appear on lengthscales which are well below our and experimental resolution (∼ 1mm), the set of parameters used in our simulations and in the experiments does not allow for a stable disorder-induced droplet, therefore sliding occurs either by precursors nucleation from the trailing edge or as first-order phase transition for top shearing. On the other hand, it is expected that in thin films (L y /L z → ∞, L x /L z → ∞), interfacial disorder may induce a droplet nucleation of the kind predicted in Ref. [25]. However, since the calculation of ∆E(r) involves two equilibrium configurations, a quasi-static model is the right candidate to tackle it.

III. Discussion
In this paper we have introduced a scalar model for the onset of frictional sliding of a three dimensional elastic object resting on a rough surface. We have devised a scalar elasticity model which allows an analytical treatment of the relevant quantities, and the straightforward implementation of the quasi-static dynamics. This model incorporates, for the first time, mesoscopic laws of contact mechanics at the frictional interface, reproducing with remarkable precision Archard's law and Amonton's first and second laws. Most importantly the scalar model is capable of reproducing with good accuracy the real contact area dynamics, the precursors' envelope dynamics prior to the transition to sliding, and the normal and shear internal stress distributions close to the slider-substrate interface. The model stems from a strong Ansatz, namely that the components of displacements u x and u z are decoupled and u y 0. However, the solution of the model is exact: if one accepts the initial Ansatz, one has at hand the analytical expression for any physical observable in static equilibrium. The numerical implementation of the model is required to take into account the statistical heterogeneity inherent to the asperity disorder of the underlying rough substrate.
The first limitation of our model has been discussed previously, and consists in neglecting the Poisson expansion and the sample torque, which can have strong implications only in the case of top uniform shearing conditions, although for very large samples the scalar models conclusions should be respected. Hence no firm general conclusions nor predictions on the occurrence of frictional sliding and precursor dynamics can be drawn based on these observations. Earthquakes faults are mostly driven uniformly from a distance implying that, in average, precursory activity should not be present. Stress gradients and hence precursor activity could, however, arise either due to local heterogeneities or because the fault plane is tilted with respect to the earth crust [29]. The lack of universal scaling forms dictating the precursors evolution, however, makes any forecasting of catastrophic events extremely difficult, especially when we do not know precisely the loading conditions. Mat. sec. F). To estimate these parameters, we only need a direct comparison with simple experiments, such as the validation of the Archard's law [30,31] to tune the normal stiffness (see Fig.2(a)), and an experiment like the one reported in Ref. [32] for the transverse stiffness.
This eliminates from the picture a host of dynamical quantities that are often difficult to quantify, or even to justify from the physical point of view. This is the case, for instance, of phenomenolgical friction coefficients interpolating between statics and dynamics employed in 1D [8] and 2D models [11][12][13], or the bulk damping coefficient γ, whose numerical value is usually put in by hand [6][7][8][9][10][11][12][13], or of the actual value of the reattachment delay time τ responsible for the frictional interface contacts rejuvenation [6,7,10]. While it is out of doubt that, upon cessation of motion, the contacts at the interface reform and strengthen [5,31], it is very hard to infer the rejuvenation characteristic time τ from experimental data. In our model, we did not include the contacts reformation process within our quasistatic protocol, showing that it is not a necessary ingredient to recover the precursors overall profile. Second, our model may assess which of the observed experimental features are due to the outof-equilibrium dynamics and which are mostly due to equilibrium properties. For instance, our model is able to recover the precursors steps and the shape of their envelope, but fails to reproduce the increase in the precursors waiting times when the shear is applied at higher and higher h. Thus, we can conclude that this intriguing aspect is probably due to inertial effects present when shear is applied through an external spring (see Eq. (22)). To check this experimentally, it would be sufficient to change the spring displacement U S rate and detect any possible change in the leaps phenomenology. To the contrary, our model allows to establish that the occurrence of precursors is in fact a quasi-static physical process. Any of the equilibrium states reached by the slider during the adiabatic evolution, is just one of the meta-stable configurations in which the system can dwell. This large number of meta-stable states is mainly due to the disorder heterogeneity of the roughness at the interface, and to a much minor degree to the rules adopted to detach the contacts when they satisfy the condition τ C (x, y, 0) > 0.
Because of its quasi-static nature, our model cannot reproduce the detachment front dynamics. According to the definition provided in Refs.[1, 2, 4, 27, 33] a detachment front indicates a drastic reduction of the real area of contact which takes place on time scales which are roughly in the millisecond range. The entire precursor experiment occurs instead over a few minutes [4,27]. Experiments have revealed three different types of crack-like rupture fronts, slow, sub-Rayleigh, and intersonic (or supershear), according to their propagation velocity through the frictional interface. Precursors advance by arrested front propagation: discrete increments, indeed, occur by rupturing the contact interface at a velocity which corresponds to sub-Rayleigh fronts at the begining, and to slow fronts close to the sliding transition [27].
In particular a final slow front is responsible of the static-to-dynamics frictional sliding. Our model does not capture the crack-like propagation of fronts, since the fronts are a dynamical out-of-equilibrium processes in between two equilibrium states, namely between precursors. Nevertheless, our model might substantiate the experimental observation on the relation between precursors appearance and slow fronts triggering the frictional sliding. As a matter of fact, in Fig.7 we were able to reproduce quantitatively the shear and normal stress profiles before any precursor nucleation occurred. In Ref. [2] these stress distributions were related to the ensuing slow rupture front (see Fig.2A in [2]). Thus it is possible to argue that whenever we observe a precursor activity, the transition to sliding is triggered by slow fronts.
To summarize the central finding of our work, three dimensional finite body scalar Green's function makes it possible to investigate the dependence of many physical observables on any sample parameter. Our results show that the evolution of the fronts depends in a non-universal way on the loading conditions and the sample dimensions and shape. Only for some loading condition, the precursors follow a curve that allows for a simple universal rescaling in terms of the sample dimension: this prediction can be experimentally checked.
Moreover we have shown that large stress gradients take place not only at the frictional interface but also within the material bulk. These gradients are mainly due to the way the external shear is applied and to the sample geometry, on top of frustrated Poisson expansion and elastic torque. Hence no firm general conclusions nor predictions on the occurrence of frictional sliding and precursor dynamics can be drawn based on these observations. Earthquakes faults are mostly driven uniformly from a distance implying that, in average, precursory activity should not be present. Stress gradients and hence precursor activity could, however, arise either due to local heterogeneities or because the fault plane is tilted with respect to the earth crust [29]. Measurements of local variations of the Coulomb stress around earthquake faults have been used to assess the correlation between stress accumulation and earthquake triggering [15][16][17]. Predicting earthquakes based on slip or stress accumulation has been so far an elusive task [18,19], and the reason behind this failure can be addressed in the scenario pictured in our analysis. Indeed, as illustrated, we find that for loading condition leading to large stress gradients, the evolution of the Coulomb stress measured above the contact interface provides only a rough indicator of the ensuing detachment front dynamics, which instead appears to be very well characterized by the real contact area variation. Furthermore, the lack of universal scaling forms dictating the precursors evolution, makes any forecasting of catastrophic events extremely difficult, especially when we do not know precisely the loading conditions.

A. The scalar model
We consider an elastic PMMA macroscopic body resting in equilibrium on a rough surface.
We derive the equations for the displacements of the elastic body subject to external forces, within the scalar elasticity approximation. In particular, we are interested in the solutions for the displacement fields at the frictional interface. There is no direct connection of the model theory to actual elasticity. The latter involves three displacement components and a system of coupled equilibrium equations enforced with boundary conditions. Let us first illustrate the physical system that our model aims at reproducing, which coincides with the experimental setup described in Refs. Since no external force is acting on theŷ direction, and because the sample geometry fulfills the condition L y L x , L z , we have assumed Moreover, the scalar elasticity yields that the stress tensor satisfies σ yk = 0 (where k = x, y, or z). Thus the scalar equations for the decoupled displacement fields take the following where E represents the Young's modulus and ν the Poisson's ratio. The scalar elasticity Eqs.
(2) can be analyically treatable, once one specifies the proper boundary conditions.
At the equilibrium, internal stresses at the surface must counterbalance the external forces 13 acting on the sample. Since we consider a slider of dimensions [L x , L y , L z ] the boundary conditions for Eqs. (2) are As shown in Fig.1(a), F lat S corresponds to a shear force in the x direction, applied to the elastic slider on a portion of the plane (0, y, z) of size L y × 2∆h centered around z = h and θ(x) stands for the Heavyside step function; F top s is a shear force (also pointing to the x direction) uniformly applied on top of the slider; F N is the normal force, i.e. a force applied on the entire top plane and pointing toward −ẑ; the surface stresses ϕ x,z surf (x, y, 0) represent the interaction between the elastic body and the rough underlying surface at the plane (x, y, 0), in the x and z direction respectively (see Fig.S1). With the boundary conditions  [34]. The result reads where G(x, ξ; y, η; z, ζ) is the Green function: with γ 000 = 0, γ 0mp = γ n0p = γ nm0 = 1/2, γ 00p = γ 0m0 = γ n00 = 1/4, and γ nmp = 1 otherwise.
In the Suppl. Mat. sec. A we furnish the analytical procedure for the fast convergence of the sum in Eq. (6).
The quantitites u x and u z represent two arbitrary constants corresponding to the displacement fields averaged over the entire sample volume. They must be chosen imposing two additive equilibrium constraints, namely that each component of the the surface forces, over the whole contact surface, must be equal and opposite to the external forces: The last step is to provide an adequate analytical expression accounting for the effects of the interface on the slider mechanics, this is done by introducing the surface forces ϕ x surf and ϕ z surf . In Suppl. Mat. sec. B we derive the expression of these forces, according to the contact mechanics theories. In first approximation they are both linear in the displacements u z and u x , i.e., for the force acting on theẑ direction, we have is a random displacement that models the height fluctuations of the rough substrate (see Fig.S1). It is formally an uncorrelated noise extracted from a Gaussian distribution with a variance σ 2 ξ = 1µm, being σ 2 ξ the roughness of the underlying surface [1-4, 31, 35] (see Suppl. Mat. sec. F). The interfacial interactions alongx are given by and u x (x, y, 0) 0. The constant c z and c x appearing in the expression of k z and k x are the only two adjustable parameters that our model encompasses (see the next subsection and Suppl. Mat. sec. F). The expressions for k x and k z respect the laws of contact mechanics [28,36,37] and are entirely motivated by experiments: the transverse (or tangential) stiffness k x of PMMA is indeed proportional to the normal load [32], which in general decreases exponentially with the vertical elastic displacement [38]; the normal stiffness is Therefore, internal stresses are not decoupled at the interface, as they are connected via the local normal pressure entering the definition of the tangential stiffness k x .
Finally, introducing the linear expressions (8) and (9) into the Eqs. (5) and (4) respectively, we obtain closed equations for the displacements at the contact plane: In the former expressions for the interfacial displacements, terms involving the contributions arising from the external shear or normal forces F S and F N can be calculated analytically. This is, indeed, one of the novelties that our model introduces: the complete expression of the Green function (see Eq. (6)) allows the determination of any of the force-induced components in the interfcial displacements equations. It will be clear in the next sections that this property entails the critical interpretation of the experimental and numerical results and, in particular, it furnishes precise predictions on the precursors' appeareance and dynamics and their dependence on the slider dimensions. In the Suppl. Mat. sec. C we give the full analytical derivation of the terms proportional to F top s , F lat s and F N appearing in Eqs. (10) and (11). Hereby, to simplify the displacements expressions, we introduce the following shorthand notations: thanks to which the Eqs. (10) and (11) take the form

B. Discretization, numerical solution and quasi-static dynamics
The slider interfacial displacements, i.e. the solutions of the Eqs. (15) and (16), are achieved by discretization of the slider bottom plane. As a matter of fact, we place a square grid on the contact plane, composed by elements having an area ∆x×∆y, so that L x = N x ∆x and L y = N y ∆y with ∆x = ∆ y = 1mm (see Fig.S1). Albeit the terms in Eqs. (15) and (16) are defined on the entire contact plane (x, y, 0), we calculate them only on each single point (x, y), which is the center of the grid element. This is the case, by instance, of the surface forces ϕ x and ϕ z in Eqs. (9) and (8) respectively, which are formally distributions: we interpret them as acting effectively only on the grid center point, representative of the enclosed area ∆x × ∆y, as shown in Fig.S1.
In the Suppl. Mat. sec. D we report the formal derivation of the discretization technique, which leads to the following expressions for the linear inversion Eqs. (15) and (16): and where the matrices A x ij and A z ij are defined in Eqs.(D19) and (D20) respectively and the vector v 0 z in Eq.(D21). With the former expressions at hand, we can calculate the equilibrium displacements alongx andẑ, compatible with given values of the external normal and shear forces: we hereby recall that the two constants u x and u z are set to ensure that the surface forces counterbalance the external loads (see Eqs.(D23) and (D25)).
In a typical simulation, external shear forces are increased quasi-statically and the actual values of the local interfacial displacements are calculated numerically at the discretized bottom interface thanks to Eqs. (18) and (17) respectively (see Fig.S1). Indeed, we first check for the equilibrium alongẑ, and secondly alongx. Contact springs are disconnected, i.e. irreversibly broken, whenever the local Coulomb stress satisfies where µ represents the static local friction coefficient set to µ = 0.5. When the condition (19) is fulfilled, the corresponding interface portion is detached from the underlying surface resulting in a local slip event. Every time a spring is broken, the equilibrium Eqs. (18) and (17)

C. Model calibration
One of the key features of our scalar quasi-static model is that Archard's principle is imposed at the mesoscopic scale [30] since the area of real contact of a grid element is function of the vertical load acting on it, i.e.
. This is in agreement with the mechanics of the interfacial asperities, whose transverse stiffnesses , and with the elastic-response picture emerging from the experiment of Berthoud and Baumberger [32]. The definition of the on-site real contact area offers the advantage of controlling the fluctuations of the total contact area during the entire quasi-static evolution. As a matter of fact, it is possible to define the the total real area of contact as and monitor its change as a function of the applied loads F N and F S . This is a notable progress provided by our model when compared with previous 1D [6-8, 10, 14] and 2D models [11,13,40], allowing for a direct quantitative comparison with experiments whose main analysis tool is, indeed, the observation of the real contact area evolution [1, 2, 4, 31].
We The normal stiffness is obtained by measuring the real contact area as a function of the normal load when no shear is applied, and tuning c z until the resulting area matches that reported in Ref. [31] (see Fig. 2(a)). Indeed, as detailed in the Suppl. Mat. sec. F, we define the total real area of contact (20) in the discrete form as Changing the value of the constant c z corresponds to change the equilibrium set of u z [i] given by Eq. (18): higher is c z , stiffer are the interfacial springs, smaller will get the coresponding real contact area. The best value for c z = 1.65 × 10 8 N m 2 yields the curve reported in Fig.  2(a), showing a remarkable agreement with the experimental observation.
To determine the transverse stiffness, we compare the quasi-static evolution of A R detected in experiments with the corresponding one obtained from simulations ( Fig. 3(a)). In particular, we consider a block of dimensions L x = 140mm L y = 6mm and L z = 75mm under a normal 19 load F N = 3.3kN and increase adiabatically the lateral force F lat S applied at height h = 6mm as in Ref. [4]. As shown in Fig. 2(b), as the lateral shear force is increased, the portion of contact area close to the trailing edge decreases drastically. According to the definition furnished in Ref. [4], precursors correspond to the regions of the frictional interface which undergo a reduction of the area of real contact, for values of the applied shear well before the static frictional force. A pictorial view of the adiabatic precursor evolution is reported in the inset of Fig. 2(b), where the color code represents the variation of the average local contact area with respect to its value at F S = 0. The boundary between the portion of contact surface which has decreased and that which has increased during the shearing process, determines the precursor size . This yieldis a curve that we compare with experiments to estimate the best value of c x = 1.65 × 10 12 N m 3 (see the caption of Fig. 2(b) and Suppl. Mat. sec. F for more details).

D. Loading mode and stick-slip events
Throughout the paper, we will consider the conceptually simple case in which the sample is loaded by imposing a constant shear force on the appropriate boundaries. This means that we will adopt F S as the adiabatic variable (quasi-static parameter), and calculate the equilibrium interfacial displacement from Eqs. (17) and (18) each time that F S is slowly increased. This leads to discrete "leapfroggy" precursors for which we study the continuum envelope as reported for instance in Fig. 2(b). However, the discrete nature of the precursors dynamics is more apparent if we drive the system as in the experiments reported in Ref. [4], where the lateral force is applied through a spring with elastic constant K S = 4 × 10 6 N/m.
To model this case, we replace the external force appearing in Eq.(4) with the expression where U S is the externally applied displacement, which now corresponds to the adiabatic adjustable parameter. u x , on the other hand, has the same meaning as the force-controlled protocol. In Fig. 4, we report the evolution of the spring force (Eq. (22)) as a function of the applied displacement for a typical simulation. Small stick-slip events, corresponding to discrete precursors leaps, are shown in the inset of Fig. 4 We can then study the different terms appearing in the former expression. The first sum is , which can be straightforwardly evaluated [42] F (x, ξ; L x ) = 1 48 The second term can be recast as [42] and we can get a compact representation introducing the function Hence Eq.(A3) acquires the final form The same can be done for the third term in the righthand side of Eq.(A1).
However a close analysis reveals that the ensuing formula does provide a fast convergence for any set of points (x, ξ; y, η; z, ζ), but for those who lie on the hypersurface x = ξ, y = η, z = ζ.

B. Surface forces
To complete the solutions in Eqs.(4), (5) for the displacements fields at the contact plane, we have to introduce a functional form for the surface stresses ϕ x,z surf . This expression represents the slider-surface interaction, embodying the microsopic details into a coarsegrained description. Albeit the surface forces ϕ x,z surf (x, y) are formally distributions defined on the contact plane (x, y, 0), they must be interpreted as effectively acting only on a single point (x, y) representative of the surface portion ∆x × ∆y that surrounds it (Fig.S1). In the next section we will provide a clear description of the slider-surface interface discretization scheme. Here we only clarify the assumptions which ϕ x,z surf (x, y) are built from: • the surface element ∆x×∆y behaves as a macroscopic object, fulfilling the macroscopic laws of friction; • the surface forces are purely elastic, i.e. linear in the displacement fields.
We consider first the surface force along the z direction ϕ z surf . Experimental evidence [38]shows that an elastic body, under a squeezing pressure P , exhibits an average surface separation u z which decreases with P as here u 0 z depends on the nature of the surface rougheness but is independent on P . Thus, the first assumption requires that the internal stress σ zz (x, y, 0) on the microscopical surface ∆x × ∆y centered around (x, y), can be written as The minus sign relates to the tensorial nature of σ zz (x, y, 0), always pointing along the −ẑ direction. Now, since at the equilibrium σ zz (x, y, 0) = −ϕ z surf (x, y), as noticed in Ref. [37] the local interfacial stiffness is given by which, for u z (x, y, 0) u 0 z (x, y) can be recast as where c z is a constant to be set. Thus, the surface stress can be written in the following linear form fulfilling the second assumption. We notice that the previous form holds only for u 0 z (x, y) > u z (x, y, 0), and ϕ z surf (x, y) = 0 when u 0 z (x, y) < u z (x, y, 0) (u z (x, y, 0) ≥ 0), see Fig.1(a). We now consider the interfacial stress along x ϕ x surf (x, y). For small displacements u x (x, y, 0), the second assumption requires that it can be expanded to first order as where u x (x, y, 0) 0. We have to find an expression for k x (x, y). From Eq.(B2) the normal load acting on the surface element ∆x∆y is given by Now, the experiment in [32] shows that an elastic body behaves like an harmonic spring (with constant K x ), if subject to a small shear force: the spring constant is linear in the applied load, i.e.
Hence, fulfilling the first hypothesis, we can assume that the surface element ∆x × ∆y is characterized by a microscopic transverse spring stiffness with c x constant to be set. Thanks to (B6) and (B9) we can finally write the interfacial stress along the x direction as Finally, introducing the linear expressions Eqs.(B10) and (B5) in Eqs. (4) and (5) respectively, we obtain the closed Eqs. (10), (11).

C. External force contributions to the interfacial displacements
In the expressions (10) and (11) for the interfacial displacements, terms involving the contributions arising from the external shear or normal forces F S and F N can be calculted analytically.
In Eq.(10) we can evaluate the term proportional to the shear force applied uniformly at the top surface F top s . According to [42], plugging Eq.(6) in Eq.(4), we obtain In the case of a shear force applied uniformly on the slider trailing edge (x = 0, h = ∆h = Lz 2 ), the former expression simplifies to Again we can define the contribution to the local displacement due to a lateral shear force We can perform the same analysis for normal displacements, obtaining again a constant contribution due to the normal force

D. Model discretization
We hereby show the procedure to discretize the contact surface between the slider and the underlying interface. Therefore we take a grid of the slider bottom plane, with an individual element having an area ∆x × ∆y, i.e.
as it is shown in Fig.S1. The central point (x, y) of each grid element takes the following discrete form The solutions Eqs.10 and 11 at any point (x, y) require the integration over the whole surface, i.e. the sum over the discrete set of the points (ξ, η) Therefore any grid center point (x, y), as well as any point (ξ, η) can be represented by a single index i and j respectively, and displacements and surface forces become to one-dimensional vectors made of N x N y elements each: Correspondingly the Green's function on the contact plane (z = ζ = 0) is a N x N y × N x N y matrix: where each element can be calculated analitically thanks to Eqs.
Finally, from Eq.(C6) we find that the contribution given to the displacement u z by the loading force F N is constant over the whole interface, so that We can now write the expressions (10) and (11) for the discrete component i of u where we made use of Eqs.(B5) and (B10). Now, let us introduce the matricesk x andk z , defined as Thus expressions (B5) and (B10) can be cast in vectorial form as 29 where we introduce the vectors u x and u z whose components are constant. Then we can obtain the solutions by inversion: whereÎ is the identity matrix. Finally we have to impose the constraints in Eq. (7) that can now be written as In order to fulfill Eq.(D18) we need the solution of Eqs.(D14) and (D15) for the i-th component of the displacement fields. Introducing the simplified expressionŝ we obtain the final expressions (17), (18). Inserting Eqs. (17) and (18) in Eq.(D18) and making use of Eq.(D12), the x constraint can be written as from which

30
In the same way we can handle the constraint along z: from Eq.(D18) we have and Let us introduce the short notations for the internal stresses at the interface After discretizing the contact plane in a mesh, i.e. τ surf (x, y) → τ surf [i] and σ surf (x, y) → σ surf [i] (see Eqs.(I2) and (I4)), the local friction law (19) requires that any site i undergo the macroscopic laws of friction. In particular whenever the discrete condition (I5) is satisfied, the corresponding spring is considered broken: Then we insert the previous relations in Eq.(D23), achieving Plugging Eqs.(E3) and (E4) in (17), Now, the requirement Eq.(E2) translates to  (D5); to the contrary, some of the formerly "attached" site could now detach: for these k x [i] = k z [i] = 0. In either case, a new equilibrium along z is needed. When the values of k x and k z in the list of available springs has not changed, we proceed further to the equilibration along x Eqs.(D23)- (17). After the equilibrium has been reached, it is still possible that one of the attached sites could now fulfill the friction law Eq.(E2) and get broken; if more than one site satisfy the rupture condition, we erase that with the largest elastic energy Hence, since the list of available springs has newly changed, we go back to the equilibration along z. It is important to notice that we do not increase further F S untill the list of available springs has not been modified by any of the possible events: attachment, detachment or rupture. When none of these occur, F S is further increased according to Eq.(E6) and the quasi-static protocol is repeated until all the springs are broken.

F. Model calibration
In this section we describe the parameters that enter our model. They can be divided into three sets: material parameters, sample parameters and adjustable parameters.
• Material parameters. To this set belong the Young's modulus E, the Poisson's ratio ν and the friction coefficient µ.
• Sample parameters. To this set belong the macroscopic dimension of the sliders L x , L y and L z , the loading force F N and the parameters concerning the lateral shear F lat S : h and ∆h.
• Adjustable parameters. To this set belong the parameters that characterize the frictional properties of the slider-bottom plane interactions in our model. They are not directly inferred from the experimental setup, since they depend on to our main assumption that surface forces are purely elastic; however we can adjust them in order to reproduce the experimental features of Ref. [31]. They are the mesh sizes ∆x and ∆y which set the degree of accuracy in the description of the contact plane slider-rough surface; c x and c z which will set the stregth of the linear interactions between the slider and the underlying surface along the x and z directions respectively; the "noise" along the z direction u 0 z , i.e. the springs rest lengths: this is independent of the slider properties but depends only on the surface roughness.
The following discussion will only concern the set of adjustable parameters. The value of ∆x and ∆y are crucial for the computing time. As a matter of fact, from Eq.(D1), the smaller the mesh sizes, the higher will be the number of elements N x and N y , increasing considerably the number of the vector and matrix components appearing in our analysis.
The values of u 0 z will be set as follows. Firstly we notice that it corresponds to the average "height" that the underlying rough surface attains within the area ∆x × ∆y centered around the point (x, y) (see Fig.S1). In general (see by instance Ref. [36]) a rough surface is selfaffine on distances ξ, i.e. given two points (x, y) and (x , y ) they are correlated whenever where ξ is the surface correlation length. As a consequence of this, the noise u 0 z (x, y) is a correlated Gaussian noise on scales ξ, but is an uncorrelated Gaussian noise for distances ξ. Throughout our analysis we take ∆x ξ = ∆y ξ 1, since we assume ξ 10µ for PMMA and having chosen ∆x = ∆y = 1mm. We can thus consider the components u 0 z [i] as uncorrelated and Gaussian distributed according to The average height u 0 z is the root mean square surface roughness Finally we consider both constants c x and c z . To set the the constant c z we refer to the experiment reported in Ref. [31]where it was shown that the real contact area increases linearly with the applied load where the constant α can be determined from the experimental curves (see Fig.1(b)), and ρ(x, y) represents the real contact area density. According to the assumption that each surface element behaves as a macroscopic object, we can assume that Eq.(F2) is valid locally Ly 0 dηG(x, ξ; y, η; 0, 0)ϕ x surf (ξ, η, 0)+ + K S (Us− ux )

Ly2∆h
Ly 0 dη h+∆h h−∆h dζG(x, 0; y, η; 0, ζ) , (G1) and adopt the shorthand notation where the value of the integral is furnished in Eq.(C3). Notice that in this case we consider only the lateral contribution to the external shearing force, neglecting, for the moment, the force exerted on top of the slider. Once we apply the discretization scheme at the frictional interface, the displacements alongx become which in vectorial form read where we implicitely made use of Eq.(D12). Inverting Eq.(G4) we finally achieve where we recall that u x is a vector with all components constant.
The force balance in Eq.D18 now takes te form By inserting Eq.(G5) into Eq.(G6) and making use of the defintion Eq.(D19) we obtain, after straightforward passages, This complete the solution, indeed after Eq.(G7) has been plugged into Eq.(G5) the set of the displacements at the frictional interfaces can be obtained.

(H9)
The second term in Eq.(H2) is straightforwardly evaluated (see [42]) Finally, the third term is In Fig.S10 we show the quasi-static evolution of the normal and shear stress at the reference plane z P = 2mm. We notice that in the limit z P → 0, Eq.(H2) and Eq.(H12) transform to E1. Finally we can define the Coulomb stress at any point in the sample τ C (x, y, z) = |τ (x, y, z)| − µ |σ(x, y, z)| . (H13) When we perform the mesh discretization at the contact plane, we have to map the variables The same can be done for the third term, although from Eq.(H11) it is apparent its dependence on x: The discrete form of the shear stress can be cast as and in the vectorial form thanks to Eq.(D12). Applying the same analysis to the normal stress Eq.(H12), we have for the second term the same as Eq.(H15), i.e.
and finally which in the vectorial form transforms to thanks to (D13). Hence the discrete expression of the Coulomb stress is given by and its quasi-static evolution is shown in Movie S2.

I. Geometrical dependence of front precursors
We report here some general consideration on the role of geometry in the shape of precursors. Plugging the continuous expression of Eq.(17) into the first of Eqs.(E1), we obtain that the shear stress at the interface can be written as or, in its discrete form, as A x (x, y; ξ, η) is the continuous version of the matrix A x ij introduced in Eq. (17), connected to the inverse of the Green's operator and to k x (x, y) (see Eq.(D19)). According to Eq.(I1), τ surf (x, y) can be decomposed into three contributions stemming from (i) forces applied to the top surface F top S , (ii) forces applied on the trailing edge F lat S and (iii) the average volume shift (see respectively Eqs.(C2), (C5) and Eq.(D23)). Among these, only the second term is responsible for the stress gradient at the interface. Indeed, when loading is exerted on the top plate, u lat S = 0 and u top S is uniform across the entire slider bottom plane as it turns out from Eqs.(C1) and (C2). Moreover u x is constant for any site, being the avergae displacement.
The normal stress at the interface, on the other hand, takes the following continuous and discrete expression respectively (v 0 (x, y)), which in average provides an uniform contribution; on the normal applied load F N , which is again even and uniform throughout the frictional interface (the edge effects are playing indeed a minor role); and on u z which is by definition uniform.
With the discrete expressions for the surface shear and normal stresses at hand, we can write the local slip condition (19) as Based on the former analysis, when shear is applied uniformly from above, i.e. F S ≡ F top S , our model suggests that both shear (τ surf ) and normal (σ surf ) component of the Coulomb stress are uniform on the contact surface and therefore no precursors are expected. This is clearly seen in Fig.3(d) and corroborates the finding of one-dimensional friction models [10]. However, scalar model does just furnish an average picture of the internal stress behaviors when shear is applied uniformly on the top, as it does not account, by instance, for the elastic Poisson expansion, nor for torque, responsibles for the stresses heterogeneities observed at the sample edges [2].
When the shear force is applied on the trailing edge, i.e. F S ≡ F lat S , only the second term proportional to u lat S survives in Eqs.(I1), (I2). Therefore the shear stress τ surf displays an evident gradient across the contact interface, leading to detachment of the regions where this gradient is more pronounced and, eventually, to precursors activity. In particular, when loading uniformly from the trailing edge (∆h = L z /2), the corresponding contribution given by u lat S to τ surf is and the interface shear stress dependence on the sample geometry is entirely encapsulated in the ratio R = Lx LyLz (see Eqs.(C4), (C5)). Since we argue that the precursors envelope curves reflect the symmetries appering in the Coulomb stress τ C (x, y, 0), the rescaled precursor profile /L x should only depend on R, in complete agreement with the numerics reported in Finally, the role of the contact plane disorder heterogeneity on the interfacial stresses appears in Eqs.(I1), (I3) through A x (x, y; ξ, η), A z (x, y; ξ, η), k x (x, y), k z (x, y) and v 0 (x, y).
For the relative narrow distribution of the heights u z 0 (x, y) characterizing the rough substrate here considered, the contribution of the interfacial disorder to the precursor nucleation is In addition, the slider bottom corners are rounded by a radius of 2 mm to avoid both stress singularities at the edges and frustrated Poisson expansion [2,31,40].
For the simulations we used the commercial FEM software COMSOL. We approximated the geometry and displacements by use of quadratic shape functions and chose the size and local refinement of the tetrahedral finite element mesh such that further refinement would give no appreciable gain in accuracy. Resulting displacements within elements were interpolated using the polynomials of the shape functions, stresses and strains were interpolated using gradients of the shape functions.

Model calibration and stress profiles
Our model encompasses the setting of only two parameters, i.e. the asperity spring stiffnesses k x and k z , as the material elastic constants have been chosen to reproduce the PMMA properties. To tune the values of k x and k z our benchmark is the experiment of Ref. [2], where the internal normal and shear stresses, σ(x, y, z) and τ (x, y, z) respectively, were calculated at the points (x, y = L y /2, z = z p ) where z p is a plane placed 2mm above the bottom surface. The slider dimensions were L x = 200mm, L y = 6mm and L z = 100mm while no shearing force was applied on the PMMA block. We have reproduced the experimental setup by means of our FEM model, and adjusted the values of k x and k z until the shear and normal stress profiles respectively were showing a satisfactory agreement: the results are displayed in Fig.S13. In particular, it is possible to see how our 3D FEM model not only furnishes an accurate quantitative agreement with the experimental curves (dashed lines), but accounts also for the high nonuniformity exhibited by normal and shear stress when no shearing is present. Indeed σ(x, L y /2, z p ) clearly displays divergences at the bottom edges due to the interface pinnning, as well as τ (x, L y /2, z p ). In this regard, in the experiments reported in [2,3], shear stress nonuniformity deriving from differential Poisson expansion frustrated at the interface, was reduced by leaving the block edges free to expand.
In our model, both normal and shear stress nonuniformities could be controlled thanks to the rounding of the block's bottom corners, providing a remarkable agreement of the stress profiles at F S = 0.
The situation considerably changes when F S = 0. Experimentally, normal stress σ(x, L y /2, z p ) shows stronger nonuniformities localized to regions near the block's edges, as the effect of the torque induced by the application of F S (see Fig.S14(a) and (c), dashed lines). This effect is much more pronounced when shear is applied at the sample's trailing edge (F S = F lat S ), by means of a 4mm-sized rod placed at height h = 6mm [2] (Fig.S14(a)). To partially avoid this byproduct, a controlled gradient in F N was applied by introducing a slight rigid tilt to the PMMA slider. This rigid tilt of the block was notably enhanced when shear was applied on the side. We could not reproduce such a tilt in our FEM model, since this loading controlled gradient leads to a highly non-linear system response, causing severe stability problems during the FEM solution. Furthermore, also the shear stress τ (x, L y /2, z p ) is considerably altered when F S = 0 (see Fig.S14(b) and (d), dashed lines).
As in the case of σ, τ nonuniformities appear more marked while shear is impressed from aside ( Fig.S14(b), dashed line) instead of uniformly on the block's top (Fig.S14(d), dashed line).
While trying to reproduce the internal stress profiles for F S = 0, our FEM model fails.
Indeed none of the simulated stresses, achieved for different values of the shearing force and different modes of shearing, exhibits a satisfactory agreement with the experimental curves, although qualitatively the observed trends are respected (see FigS14). We believe that our model failure is ultimately ascribable to the impossibility of transposing numerically the tilting block expedient employed in the experiment. As a matter of fact, when F S is applied uniformly on top of the slider, the agreement between experimental and numerical stress profiles appear to be slightly better than when shear is imposed from the edge.   to the experimental setup (Ref. [2]). Stress calculation was performed at a value of the shear force F S right before the nucleation of the precursor (see Fig.S10 and Fig.8(a)). L x = 200mm, L y = 7mm, L z = 100mm, F N = 6.25kN. (b) Coulomb stress calculated on the reference plane above the contact interface and averaged over the y direction (solid green line): τ C (x, z P ) = |σ(x, z P )| − µ |τ (x, z P )|. Comparison with data inferred from Ref. [2] is excellent (green symbolsdashed line). The magenta line represents the y-averaged Coulomb stress at the frictional interface (τ C (x, 0) = Ly 0 dy τ (x, y, 0)): although no detachment front is yet present at the contact plane (τ C (x, 0) < 0 throughout the surface), the value of the Coulomb stress at the reference plane z P = 2mm can exceed locally the threshold, leading to the erroneous conclusion that Amonton-Coulomb law might be violated.  Eq.H12 (Eq.H20 in its discrete form) and averaged along the y direction, i.e. |σ(x, z P )| during a shearing protocol which involves both an edge and top pulling. The dashed black line corresponds to the F S value prior to precursor nucleation: the corresponding shape of |σ(x, z P )| is shown in Fig.7(a) (red curve). (b) Internal average shear stress obtained from Eqs.H2-H17 and averaged along the y direction, i.e. |τ (x, z P )|. The dashed black line corresponds to the F S value prior to precursor nucleation: the corresponding shape of |τ (x, z P )| is shown in Fig.7(a) (blue curve).    and experiments is not satisfactory (although slightly better than lateral shearing, in panels (a) and (b)). Dotted lines represent the normal and shear stress profiles calculated according to the scalar model. In this case the stresses appear uniform as a consequence of the decoupling of the displacements u x and u z (scalar elasticity approximaton).