A parameter-free mechanistic model of the adhesive wear process of rough surfaces in sliding contact

In order to develop predictive wear laws, relevant material parameters and their influence on the wear rate need to be identified. Despite decades of research, there is no agreement on the mathematical form of wear equations and even the simplest models, such as Archard's, contain unpredictable fit parameters. Here, we propose a simple model for adhesive wear in dry sliding conditions that contains no fit parameters and is only based on material properties and surface parameters. The model connects elastoplastic contact solutions with the insight that volume detachment from sliding surfaces occurs in the form of wear particles, the minimum size of which can be estimated. A novelty of the model is the explicit tracking of the sliding process, which allows us to meaningfully connect particle emission rates and sizes to the macroscopic wear rate. The results are qualitatively promising, but we identify the necessity for more controlled wear experiments and the parameters needed from such work in order to fully verify and improve our model.


Introduction
In engineering practice, wear is usually predicted by (semi-)empirical laws instead of mechanistically informed models. Unfortunately, these laws are numerous and not necessarily transferable, while experiments which reveal sufficient information to validate them are scarce (Meng and Ludema, 1995). In order to further our fundamental understanding of the wear process, parameter-free wear predictions based on microscopic mechanisms are needed. The main obstacle in this path is the multitude of those mechanisms (reviewed in Vakis et al., 2018), spanning from corrosion, local melting, subsurface plasticity, and phase transformation to fatigue and fracture from the nanoscale to the macroscale. A universal, multi-physics description of this process is thus currently out of reach and scientific progress in this area relies on describing idealized wear modes.
Hence, our focus lies here on the case of adhesive wear of brittle materials in the unlubricated regime. Adhesive wear is perhaps one of the most common types of wear, in which wear particles detach at points where strong enough junctions form between asperities from two contacting surfaces sliding relative to one another (Rabinowicz, 1995). The earliest works on wear already recognized that the wear volume V -i.e., the loss of volume of relatively sliding surfaces-is proportional to the applied normal force F N and the relative sliding distance s, at least in a certain load range (Reye, 1860;Rabinowicz and Tabor, 1951). The value of interest is therefore the prediction of the wear rate W = dV /ds. Bowden and Tabor (1939) provided the critical insight that all surfaces are rough and that the value of the real contact area thus depends on F N . With this, Holm (1967), Burwell and Strang (1952), and Archard (1953) developed wear equations of the form which assume that the wear rate is proportional to the real contact area given by F N divided by the material's hardnessp with some proportionality constant k. The contribution of Archard (1953) is the insight that wear progresses due to the detachment of wear particles at contact junctions between the surfaces instead of by atom-by-atom attrition. Nevertheless, the predictive power of Eq. 1 is low, since the constant k, commonly called "wear coefficient", remains an empirical parameter whose values can span several orders of magnitude. While Archard (1953) treated the wear particle formation as a random event with k as the detachment probability, recent computer modeling efforts (Aghababaei et al., 2016(Aghababaei et al., , 2017Aghababaei, 2019;Brink and Molinari, 2019) confirmed and expanded on an old hypothesis by Rabinowicz (1958) that there is a critical length scale for wear particle formation. If contact junctions are smaller than this critical size, they do not detach from the surface, but show some sort of plastic activity. The model is based on the amount of energy that can be stored locally in a contact and is thus available to propagate the crack that leads to particle detachment. The critical length scale relates to a value for the junction diameter d and can be expressed, based on the material's properties, as (Aghababaei et al., 2016) d where f is a geometry factor (often set to ≈ 3, see Aghababaei et al. (2016) and Brink and Molinari (2019) for a detailed discussion and derivation), G is the fracture energy, andτ 2 /2µ is the elastic energy density of an asperity with shear modulus µ when it is loaded to its shear limitτ . It should be noted that this concept has only been tested with relatively brittle material models and an extension to materials with ductile fracture behavior is outstanding, hence our limitation to the brittle regime in the present paper. A first approach to advance Archard's model based on this deterministic particle formation criterion was developed by Frérot et al. (2018) using elastic contact solutions of synthesized, self-affine, rough surfaces. Here, k is no longer a detachment probability, but related to the fraction of asperities with d ≥ d * . This model is revisited in Sec. 2, so we simply summarize that it does not reproduce a wear coefficient that is constant over a large load range as experimentally observed (e.g., by Archard and Hirst, 1956). Another approach by Popov and Pohrt (2018) is also based on elastic contact solutions of synthesized, self-affine, rough surfaces, but tries to avoid defining an asperity or contact junction directly, instead applying the energetic criterion underlying Eq. 2 to arbitrary areas of the surface. This model predicts that infinitely large wear particles can form and yields superlinear wear rates even when imposing an upper limit on particle size. The reason is (most likely) that having sufficient elastic energy available is only a necessary condition for fracture and wear particle detachment; as already implicit in the fracture model of Griffith (1921), there must also always be a stress concentration to propagate a crack.
Based on these results, it is not yet clear if the critical length scale of Eq. 2 is helpful in predicting macroscopic wear rates. Indeed, the value of d * decreases with increasing hardness of the material, which might suggest that a harder material can more easily produce wear debris even from small contacts and thus wear more, in contradiction to the typical experimental observation. The mechanistic model by Frérot et al. (2018) shares this shortcoming: Wear rates increase with decreasing d * .
The goal of the present contribution is therefore to show that a critical ingredient of the wear process was not considered in sufficient detail before: The surfaces slide relative to each other during wear. The explicit inclusion of a very simple facsimile of the sliding process into our model allows us to make sense of the seeming contradiction between the microscopic length scale d * and the macroscopic wear laws. Furthermore, we find that the material properties used to calculate d * cannot be divorced from the physics of the contact solution and that the relation between wear particle formation and wear rate is nontrivial. A comparison with an experimental work by Kim et al. (1986) that contains sufficient data to feed our model yields some qualitative agreement, but ultimately more detailed experiments are needed to develop predictive wear models.

Complementing a static model with plasticity
Our proposed new wear modeling approach is based on the work of Frérot et al. (2018), the basic ingredient of which is an elastic contact solution. Figure 1(a) illustrates the steps to obtain such a solution: First, two rough, fractal surfaces are generated for each simulation, representing the two sliding bodies. A fractal description is chosen because self-affinity is often observed in experiments of rubbing surfaces (Majumdar and Tien, 1990;Persson et al., 2005;Persson, 2014;Davidesko et al., 2014;Candela and Brodsky, 2016), as well as in simulations (Milanese et al., 2019(Milanese et al., , 2020Hinkle et al., 2020). The surfaces are square and periodic in x and y direction and represent half spaces in z direction. Periodic synthetic surfaces are generated in the Fourier domain with a fixed power spectrum and a uniform distribution of phases (Wu, 2000). Since all efficient contact algorithms describe a rough-on-flat geometry, the model uses Johnson's assumption that the rough-on-flat solution is equivalent to the rough-on-rough solution, provided that the former is solved using the gap of the latter as the rough surface (Johnson et al., 1985). Now, for a given contact solution at normal load F N , all contact junctions are identified using a flood-fill algorithm. Assuming that the contacts are approximately circular, the critical length scale is translated into a critical area A * = πd * 2 /4 and junctions are divided into plastically deforming and wear particle forming based on their size (see Fig. 1(b)). The idea is now that Archard's wear coefficient indicates how many junctions can form particles. The wear rate for each junction is obtained using Archard's estimation that a junction fully detaches after a sliding distance roughly equal to its diameter (Archard, 1953), therefore obtaining a wear rate equal to volume divided by diameter, i.e., proportional to the junction's contact area. Note that this means that the sliding distance needed to form a wear particle is different for each junction. Archard's wear coefficient is thus calculated as where A i is the area of contact junction i and the sum is over all junctions. The wear rate is then simply W = kω i A i , where ω is an average shape factor on the order of unity that accounts for the local contact and wear particle geometries. We call this model "static" because the relative sliding motion is only implied in the derivation of the per-junction wear rate. It is clear that a purely elastic contact solution is not realistic (Bowden and Tabor, 1939;Weber et al., 2018) and our first step to improve upon the above model is to introduce an elastoplastic contact solution. While at least full von Mises plasticity is desirable for realistic contact calculations, the required computing time is prohibitive for a large number of simulations (Pei et al., 2005;Frérot et al., 2019c). In order to keep the computational demands reasonable, we thus opted for a saturated-pressure solution as described by Almqvist et al. (2007) and as implemented in tamaas (Frérot et al., 2019b). The model relies on the common assumption that in elastic-plastic contact the surface pressure should not exceed the hardness of the material. To this end, saturated-plasticity approaches used in rough contact (e.g., Akchurin et al., 2016;Weber et al., 2018) solve an elastic contact problem with the additional constraint that the contact pressure p should be smaller or equal to p m , the flow pressure. Often, this flow pressure is taken to be equal to the indentation hardnessp, which can in turn be estimated from the material's uniaxial yield strengthσ asp ≈ 3σ under the assumption of spherical indentation . Note that this is an approximation of the full elastoplastic contact solution and the readers are referred to Frérot et al. (2019a) for a comparison between von Mises plasticity and saturation pressure models. We start our model development by investigating the influence of plasticity on the wear prediction for the static wear model. For this set of simulations, we describe all quantities in reduced units in terms of the Young's modulus E and the side length L of the surface. We used five realizations of each pair of rough surfaces, each surface being discretized into 4096 × 4096 pixels and described by a Hurst exponent H = 0.7 (a typical value observed experimentally, see Persson et al. (2005)), a lower wavelength cutoff of λ l = L/256, an upper wavelength cutoff of λ u = L/16, and a root mean square (RMS) of heights of h RMS = 0.005 L. For  The models in this work are based on the elastoplastic contact solution of two rough surfaces, which is approximated by solving the problem of the contact between a rigid surface (whose geometry is equal to the gap between the two original surfaces) and a flat elastic-perfectly plastic half space (Johnson et al., 1985). From this, we obtain the contact junctions shown in yellow at the bottom. (b) A first model proposed by Frérot et al. (2018) states that Archard's wear coefficient k is equal to the ratio of total area of the junctions large enough to form wear particles divided by the real contact area. (c) Our improved model takes the sliding process into account. The top surface is shifted by a small length δs each step. Similar to Frérot et al. (2018), we mark junctions that form wear particles. In contrast to the earlier work, we ignore junctions that have formed a wear particle in a previous step. All others contribute a total wear volume δV each step. After a running-in period, the wear rate can be estimated as W = dV /ds. Note that the δs between the two contact maps shown here is quite large in order to make the evolution more easily visible. each of those pairs, 1280 different offsets parallel to the surface planes between the top and bottom surface were used, since they were needed anyway for the sliding model described in section 3 later on.
We used an isotropic material with Poisson's ratio ν = 0.3. For the shear strength, values ofτ ∈ {0.006, 0.008, 0.010, 0.012} E were chosen, from which the indentation hardness was derived as p m =p = 3 √ 3τ . Values for the fracture energy were G ∈ {0.25, 0.5, 0.75, 1.0} · 10 −6 EL. The geometry factor was set to f = 3 under the assumption of spherical asperities, leading to a range of d * /L ≈ 0.004 to 0.064, i.e., in between λ l  Frérot et al. (2018) with the addition of a saturated-pressure contact solution instead of an elastic one. The wear rates are normalized by the shape factor ω. The influence of varying (a) toughness and (b) hardness is shown. The colored area represents the standard deviation of the data between different simulation runs. and λ u . Sixteen different normal loads were used, so that in total 409,600 contact solutions were computed (the fracture energy only appears in d * and has no bearing on the contact solution).
Despite the addition of more realistic contact physics, the results resemble the work of Frérot et al. (2018), see Fig. 2. The wear rate is superlinear and the wear coefficient increases with load, tending towards a limit of one. Additionally, it should be noted that while the wear rate increases with decreasing fracture energy (first row of Fig. 2), which seems reasonable, it also increases with increasing hardness (second row of Fig. 2), which contradicts Archard's wear law. This is the same in the original work and is a simple consequence of using d * as a wear particle formation criterion: The bigger d * , the fewer junctions on a given contact map form wear particles. Since d * ∝τ −2 , this effect cannot be offset by the increasing real contact due to plasticity, which is less than quadratic (Bowden and Tabor, 1939).
The intermediate conclusion is either that the critical length scale model contradicts experimental evidence or that a wear model based on a single, static contact solution is missing a critical ingredient. We will argue for the latter by implementing a simple sliding process.  Figure 3: Running-in and fitting the wear rate. The wear rate is the slope obtained by linear regression. The initial points indicated in the inset corresponds to the instantaneous wear volume V 0 , which is much higher than any subsequent increment δV .

Sliding Model
In order to explicitly include the sliding process in the model, we extend the model presented in the previous section by continuously displacing the top and bottom surfaces against each other. We do not calculate wear rates directly from the real contact areas, but instead estimate wear volumes under the assumption that the junctions are circular and the resulting wear particles spherical (take A i as the area of the contact spot i and use d i ≈ 2 A i /π, V i ≈ πd 3 i /6). The algorithm-also shown in Fig. 1(c)-works as follows: 1. Solve the contact, find all contact junctions, and record the instantaneous wear volume V 0 as 2. Displace the top surface by a sliding increment δs and recompute the contact solution. Due to the periodic boundary conditions, the part of the top surface that leaves the simulation box on one side is in fact reinserted on the opposite side. 3. Find all contact junctions. Iterate over them and mark any junction that shares at least one pixel with a junction that has formed a wear particle in the previous step as "disabled". The same for any junction that shares at least one pixel with a junction that was disabled in the previous step. All other junctions are marked "active". This mimics the ongoing detachment of a wear particle, whose volume should only be counted once, but which still contributes to carrying normal load. Because the long-distance rolling of a wear particle cannot be modeled in the present framework, we assume the particle is lost quickly. 4. Increment the wear volume by 5. Unless the desired total sliding distance was reached, go to step 2.
The output of this algorithm is a cumulative wear volume as a function of sliding distance. We chose a total sliding distance of L/2. 6 Figure 3 shows that the wear volume increases linearly with the sliding distance after a running-in period, in which the increase can be steeper. As in virtually all experiments (Reye, 1860;Rabinowicz and Tabor, 1951;Burwell and Strang, 1952;Archard, 1953;Archard and Hirst, 1956), our wear rate prediction (obtained through linear regression as W = dV /ds) reaches a constant value, yielding a first confirmation that the sliding implementation is reasonable. Note also that the instantaneous wear volume V 0 is much higher than the increments δV . We expect this effect to be very hard to observe experimentally, since it takes place at the exact onset of sliding. Nevertheless, it has an analogy in the well-established observation that friction at rest is increased compared to friction during sliding. Moreover, running-in periods with higher wear volume are well known from experiment (Queener et al., 1965), although they can have many causes not included in our model, such as changes of the surface roughness parameters.
For the following study, we defined the sliding distance after which the system is definitely run in as L/4, which we confirmed by inspecting the simulation data. We note that the length of the running-in period should be dependent on the surface geometry and parameters, but here this length was simply chosen as a safety margin without physical insight. We chose a sliding increment δs corresponding to 4 pixels in the running-in period and 1 pixel afterwards. 3 Video 1 shows an example simulation of a smaller surface. The surface size was reduced to make the features of the model more easily visible.
Additionally, we verified that the ratio λ l /∆x, where ∆x is the side length of a pixel, gives an acceptable discretization error. The details can be found in Appendix A.

Role of the critical length scale
Some properties of the model are immediately apparent: Contact junctions wear off as soon as they reach a critical size and thus the wear volume increment δV is made up of wear particles of roughly size V * ≈ πd * 3 /6. This reveals the role of the critical length scale as not only the minimum size of wear particles, but probably also the common size. Thus infinite particle sizes, as featured in the work of Popov and Pohrt (2018), do not occur in reality because contact junctions simply wear off before they can grow that big. It is known from experiments, e.g., by Rabinowicz (1953Rabinowicz ( , 1995 or more recently by Kirk et al. (2019), that there is a distribution of particle sizes. This can have two reasons. On one hand, d * is not uniquely defined on a real surface due to fluctuations of the adhesion and local geometry, which are both expected to affect d * . On the other hand, once particles have detached, they continue to evolve and grow, losing the connection to their initial size (Milanese et al., 2019;Aghababaei, 2019). Given that the initial particle size is large compared to its growth rate (Milanese et al., 2019), we neglect any wear due to particle growth. The wear rate can thus be defined by both d * and the rate at which new contact junctions appear and reach d ≥ d * during sliding: Here, ω (= π/6 in the present work) is again a shape factor accounting for the wear particle shape and N * + is the total number of wear particles formed during the sliding distance s, which is also the cumulative number of junctions that have newly grown to a size of at least d * . Note that the total number of junctions greater than d * , denoted N * , is roughly constant after running in. For more details, see Appendix B.

Apparent contact area
Typically, it is assumed that the wear rate, just like the real contact area, is independent of the apparent contact area (Archard, 1953). We tested if this is the case in our model (see Appendix C for the details) and found for two surfaces with different side lengths that the wear rates are comparable up to a given load, after which the smaller surface exhibits a breakdown of the wear rate with increasing load. This is not related to a breakdown of the total number of contact junctions-which strongly depends on the surface size-but to a breakdown of the scaling of N * with F N , i.e., the number of contact junctions greater than d * . This seems to have direct implication on dN * + /ds and thus the wear rate. Here, we treat this as an effect of the limited size of our surfaces and regard the data after the "breakdown" as unreliable. If such a limit to Archard-like scaling with F N is a real effect remains to be investigated.

Influence of the lower wavelength cutoff and the RMS of heights
In elastic contact solutions, the real contact area scales linearly with the inverse of the RMS of slopes (Bush et al., 1975;Hyun et al., 2004). Because the latter is strongly dependent on λ l , the true contact area is sensitive to local fluctuations at short wavelengths and in extenso to the resolution down to which the system is measured or modeled (Ciavarella et al., 2000;Persson, 2001;Ciavarella and Papangelo, 2017). This is a problem for elastic contact simulations, since it becomes necessary to include the shortest wavelength of roughness that exists in the system, which might be much smaller than d * and which is likely unknown due to experimental measurement limitations . However, in elastoplastic contact solutions, plasticity can "smooth out" small wavelengths (Pei et al., 2005). We found that this also is the case in our model, see Appendix D. While it is hard to generalize to other elastoplastic models, it at least means that the present model is robust against (in practice) immeasurable small disturbances of the rough surfaces.
Another surprising property is that the wear rates do not change between two values of the RMS of heights h RMS tested by us (see Appendix E). For an elastic contact solution this would be unexpected, since the RMS of heights also influences the RMS of slopes and thus the contact solution. In our case, it seems that the nature of the sliding model combined with the saturated-pressure contact solution cancels any influence of h RMS , at least in a certain range of h RMS and F N . This can be rationalized by first taking into account that sharper peaks also increase the local pressure and thereby provide the necessary energy to flatten them, as well as the fact that the wear rates in our model depend more strongly on the number of contact junctions, which is likely not affected as strongly by h RMS . Further study is necessary to draw definite conclusions.

Reduced units
In a first step, we used the surfaces described in section 2 as input to our sliding model with identical variations of material parameters. The results are shown in Fig. 4 and Supplementary Figs. S.1 and S.2. The first observation is that wear rates and wear coefficients are an order of magnitude smaller than in the static model (Fig. 2). They reside in the high end of reported experimental values for dry sliding experiments of ceramics (Zum Gahr, 1989;Hokkirigawa, 1991;Wang and Hsu, 1996) or of steel on steel (Archard and Hirst, 1956), which range from 10 −6 to 10 −1 . Unlike the results shown in section 2, the wear coefficient here does not asymptotically tend to unity.
First, we keep the material hardness constant and vary the fracture energy G (Fig. 4(a)). This means that for the same contact solution (which does not depend on G), d * ∝ G. With larger critical length scale, the wear rate is higher when it reaches a roughly linear slope at higher loads, but the onset of wear is delayed as can be seen by observing the change of the wear coefficient k at low loads.
Second, by keeping the fracture energy constant and varying the material strengthτ , we can observe the effects of simultaneously changing d * ∝τ −2 and the contact solution ( Fig. 4(b)). Again, larger d * leads to higher wear rates, but a delayed onset of wear, this delay being weaker than in the first case. The reason is that the reduced strength implies a reduced indentation hardness and therefore more plastic activity in the contact junctions, which will thus be larger even at lower loads.
Finally, by keeping d * approximately constant and varyingτ and G together (Fig. 4(c)), we see that the wear coefficients ultimately tend towards the same value, while the onset of wear can be somewhat delayed by making the material stronger and tougher.
It is also instructive to look at the number of wear particles, or rather the rate at which they are emitted over a given sliding distance: dN * + /ds. This is shown in Fig. 5 for the same set of simulations. It should be noted that while the hardest material and the least tough material wear the least, they emit the most wear particles. Since these particles are quite small in these cases, though, the wear rates are also small. This means that more severe wear is due to the detachment of large wear particles, which is also sometimes observed experimentally (Zhang and Alpas, 1997). Of course, when d * is roughly constant (Fig. 5(c)), the wear particle formation rate and the wear rate exhibit the same trends. Thus, if one wants to minimize the emission of fine particles, our current model suggests to use a soft but tough material.
That experiments often exhibit wear coefficients at least one or two orders of magnitude below the present results (Archard and Hirst, 1956;Zum Gahr, 1989;Hokkirigawa, 1991;Wang and Hsu, 1996) can likely be explained by the reduced adhesion even in dry conditions. We assumed that the junctions have bulk strength, although surface contamination is always present and will reduce the adhesion between two bodies in contact. The effect would be to increase d * (Brink and Molinari, 2019; Aghababaei, 2019), while not affecting the bulk hardness and contact solution. This is equivalent to increasing G while keeping all other material parameters the same (see Eq. 2). Figure 4(a) shows that this pushes the onset of wear to higher normal loads, since junctions with size d > d * are less numerous at low F N . By significantly reducing adhesion, this onset can be pushed to very high loads, see for example Supplementary Fig. S.1(a). As noted in Brink and Molinari (2019) and Milanese et al. (2020), though, the critical junction size becomes geometry dependent in the reduced adhesion case and some wear particles can be formed, thereby resulting in a reduced, but non-zero wear coefficient. Unfortunately, the actual microscopic junction shear strength for reduced adhesion is not easily accessible experimentally and the above hypothesis remains to be tested against realistic data. Kim et al. (1986) A direct comparison to experimental data is not easy, since we could not find any work which simultaneously reports the surface roughness parameters, material properties, and wear rates for dry sliding conditions. Nevertheless, we could use the data for Si 3 N 4 , SiC, and Al 2 O 3 from Kim et al. (1986), who at least report the RMS of heights for their surfaces in addition to material properties and wear rates (see Table 1), even if only for rolling friction experiments.

Using data from
In terms of the surface geometries, we used a side length of L K = 200 µm (roughly the apparent contact area of the experiment, although their contact area is an ellipse due to Hertzian contact); H = 0.7, λ l = L K /512 = 0.39 µm, and λ u = L K /8 = 25 µm. As such, d * lies in between λ l and λ u , keeping the computational demands (discretization) reasonable. The RMS of heights was chosen to be the steady-state  Kim et al. (1986). All synthesized model surfaces had a side length of L K = 200 µm, Hurst exponent of H = 0.7, λ l = 0.39 µm, and λu = 25 µm. The reduced units lie in the same range as before.   Hokkirigawa (1991). The large absolute differences between simulation and experiment are thus expected due to the different setup.
value of the wear experiment (see Table 1). Note that the values forτ /E and G/EL K lie in the same range as the values we used for our generic simulations with reduced units, verifying that the choice of parameters was reasonable, at least for ceramics. The results are plotted in Figs. 6(a)-(b) and the wear coefficients are in the same range as the previous simulations (Fig. 4)-as expected. The authors of the experimental work proposed a wear law based on "contact severity", which does not seem to yield insights beyond Archard's wear law. It is discussed for completeness's sake in Appendix F. Apart from this, Fig. 6(c) shows that our wear coefficients are at least two to three orders of magnitude higher than the experimental ones (k exp = 6 · 10 −7 to 1.5 · 10 −3 ), although this is expected since the experiment was a rolling friction setup, which often has such a large difference in k compared to sliding friction setups (Hokkirigawa, 1991). The hierarchy of wear resistance matches the experimental observation: Al 2 O 3 wears the most and Si 3 N 4 the least for the presented load range. While the crossover between Al 2 O 3 and SiC at low loads observed in the simulation is also somewhat suggested by the experimental trends, the available measurement data is insufficient to affirm a correspondence. In contrast, the wear rate of Si 3 N 4 predicted by our model increases towards the higher loads in accordance with the results of section 4.1, since it has the highest d * . It is unclear if that behavior of our model is realistic or not, since we cannot translate the load ranges in sliding to the equivalent load ranges in rolling friction. It is quite possible that the typical experimental load ranges never reach the point where the wear rate of Si 3 N 4 picks up and surpasses the other materials. Some hints can be found in the sliding wear experiments on different ceramics of Adachi et al. (1997). They report wear rates up to the order of 10 µm 3 /µm in dry conditions for normal loads in the range of 2 N to 100 N, which corresponds roughly to our results. Instead of normal loads, data is presented as a function of "contact severity" (a compound value of different, independent parameters), which makes further comparison impossible.
Ultimately, we hope to inspire more detailed experiments that adequately report on material parameters, surface roughness, and wear rates over large load ranges in order to verify and improve our model. Recent advances in the mass resolution and wear mapping of sliding experiments, as for instance proposed by Garabedian et al. (2019), can be helpful in this context because they allow observations at the onset of wear and can resolve asperity-level mechanisms.

Outlook on wear modeling
Due to the computational demands, the present model is relatively simple. It is clear that a saturatedpressure contact solution is not as satisfactory as a more realistic von Mises plasticity (Frérot et al., 2019a), but existing approaches (Pei et al., 2005;Frérot et al., 2019c) are computationally too demanding for the large number of contact solutions required at present. Similarly, the spectrum of the fractal roughness is not very broad, once again to keep computational demands low. Finally, Johnson's assumption is likely invalid for elastoplastic contact and true rough-on-rough contact algorithms should be developed. Mechanistic wear modeling could greatly benefit from any advances in these areas, if only to verify the validity of the approximations used here.
There are also several physical processes that are not explicitly included. First, a reduced adhesion between the surfaces, either due to surface contamination or lubrication, is known to modify d * . However, a purely geometric treatment of d * as used here would additionally need the contact angles , which are not well defined when using Johnson's assumption. Perhaps it would be more fruitful to work on including the fracture process directly in the simulation instead, as attempted by Frérot et al. (2019a). Moreover, coupling of a fracture model to the contact solution could inherently include asperity interactions Pham-Ba et al., 2020), which are expected to greatly affect the wear rates at high loads. Finally, we assumed that the debris particles are somehow evacuated quickly and no extended third body forms, which is of course an oversimplification in many cases.
Despite all of this, a major advantage of the current work is the simplicity of Eq. 4. The main complication is to provide a reliable estimate of N * + or dN * + /ds. If a closed form solution could be obtained, the computational demand would be significantly reduced and the model could be widely and cheaply tested and applied.

Conclusion
We developed a wear model based on a combination of a simple elastoplastic contact solution algorithm, an evolution of the contact interface through a sliding process, and a critical length scale criterion for wear particle formation. The model contains no fit parameters and yields at least qualitatively reasonable results. We find that the critical length scale d * , i.e., the minimum size a contact needs to have to detach and form a wear particle, is a necessary ingredient in this model: It defines not only a minimum wear particle size, but rather a typical wear particle size. Often, this leads to materials with larger d * wearing more, which are materials that are soft but tough (see Eq. 2). Nevertheless, for smaller loads, increased toughness can suppress wear particle formation since no large enough contact junction forms. The statistical distribution of contact junctions that form is thus also a relevant parameter. It remains to be seen if the behaviors at different load ranges can be related to experiments.
Nevertheless, an Archard-like wear behavior is approximately recovered and the qualitative agreement with an experimental work is promising. However, the present comparison is not quantitatively satisfying, since we can only compare to rolling wear, while we intrinsically model sliding wear, for which we could not find sufficient data. It is our hope that such modeling inspires more detailed dry, adhesive wear experiments, which report simultaneously the necessary material parameters, surface roughness description after running in, and resulting wear rates. Preferably the work would be performed over large load ranges in order to evaluate and improve models such as the one presented in the present paper.

Acknowledgments
This work was supported by École polytechnique fédérale de Lausanne (EPFL) through the use of the facilities of its Scientific IT and Application Support Center.

Appendix A. Discretization
In order to verify that the discretization errors are sufficiently small, we ran simulations with a ratio λ l /∆x = 16 as used in the results in this paper, as well as λ l /∆x = 32, i.e., a two times finer discretization. All test simulations, presented in Fig. A.7, used roughness parameters of H = 0.7, λ l = L/256, λ u = L/16, and h RMS = 0.005 L, as well as material properties of ν = 0.3,τ = 0.010 E and varying G as indicated in the figure. We found no difference between the two settings within the error of the simulation and thus opted for a value of λ l /∆x = 16 in the interest of reducing computational demand due to discretization.

Appendix B. Number of wear particles
Equation 4 indicates the importance of the number of wear particles N * + that are created over a given sliding distance s. Figure B.8 presents a detailed look at this data for one the main simulations with τ = 0.006 E and G = 0.25 · 10 −6 EL. In Fig. B.8(a), we see the number of wear particles formed N * + (s, d * , F N ), the cumulative number of junctions N * − (s, d * , F N ) that have shrunk below d * during a sliding distance s, and the number of junctions N * (s, d * , F N ) greater than d * at a given sliding distance s. It is N * = N * + − N * − , wherefore N * − closely follows N * + after some steps. Just like the wear rate, the wear particle formation rate is linear with sliding distance. The slope dN * + /ds is obtained by linear regression, cf. the wear rate in Fig. 3. Figure B.8(b) demonstrates that Eq. 4 with ω = π/6 slightly underestimates our simulation results obtained by directly summing the volumes of wear particles, even if the simulation also assumes spherical particles. The reason is that-due to discretization-contacts may grow a little bit bigger than d * before they detach. In reality, particles would not be perfectly spherical and adjustments to ω have to be made in any case.

Appendix C. Apparent contact area
In order to study the behavior of our model as a function of the apparent contact area, we also produced five pairs of surfaces with a side length of L 1 = 0.5 L (we will not normalize material properties and results to L 1 , but stay with L as a length unit for easier comparison). Otherwise, all surface and material parameters were the same as for the surfaces with side length L: H = 0.7, λ l = L/256 = L 1 /128, λ u = L/16 = L 1 /8, and h RMS = 0.005 L, as well as material properties of ν = 0.3,τ = 0.010 E, and G = 0.50 · 10 −6 EL. Figure C.9(a) shows that the wear rates are indeed comparable up to normal loads of F N ≈ 0.5 · 10 −3 EL 2 . At higher loads, the larger surface continues a roughly linear increase of wear rate with load, while the wear rate of the smaller surface breaks down. Since the wear rate depends on the number of junctions growing to the critical size, we also looked at the average number of contact junctions as a function of load (Fig. C.9(b)). This analysis proves to be not very enlightening on its own, though: The larger surface has more contact spots, but does not wear more, and there is no crossover in behavior as a function of normal load. A comparison of the number of contact junctions N * that are larger than the critical size reveals more. In Fig. C.9(c) we can see that N * breaks down from its approximately linear scaling at the same time as the wear rate does. We interpret this to mean that smaller surfaces statistically have less contact spots overall and thus cannot accommodate the necessary number of large contact spots to keep up the linear growth of wear rate with load. Video 2 shows a simulation in the high-load regime of an even smaller surface and indicates that contact junctions also start to merge, which additionally reduces the number of individual contacts. This might well be connected to a transition to severe wear, since the contacts are expected to interact if they come closer together, at which point their wear particle sizes are no longer predictable from the individual contact areas Pham-Ba et al., 2020). Such phenomena are not included in the present model, and we have to treat data after this "breakdown" as unreliable.
Note that we chose the larger surface to obtain the results in the main text since it can sustain a larger range of loads. Because we originally started our investigations with the smaller surface, though, it was used for the verification studies in the appendices unless otherwise noted.

Appendix D. Lower wavelength cutoff of the surface spectrum
We studied the influence of the lower wavelength cutoff of the surface spectrum ( Appendix F. Wear law based on "contact severity" For the sake of completeness, we will shortly comment on a proposed wear parameter S c in the publication of Kim et al. (1986) In contrast to Archard (1953), the hardness does not play a role, but the toughness does. The experimental data in Kim et al. (1986) consists of only two loads per material and Eq. F.1 is fitted to a log-log plot, where it works approximately, with values of α = 1.56 × 10 −5 mm 3 /mm and n = 5.46. The same is true for our simulations, see Fig. F.12, although with fit parameters of α = 1.233 × 10 −4 mm 3 /mm and n = 0.9 and some deviation for Si 3 N 4 . It is not clear if Eqs. F.1 and F.2 are predictive, especially since we do not find a strong dependence on h RMS , which seems to have been included for the convenience of having a unitless S c .
10 −2 10 −1 10 0  Figure F.12: Relation between our simulated wear rates and the proposed "contact severity" parameter Sc (Kim et al., 1986). The dashed line represents the best fit of all data.
The influence of hardness is completely ignored. Indeed, later publications did not find a clear correlation between S c and W and rather propose that the maximum stress p N leads to the correlation, similar to F N in Archard's wear law (Hsu and Shen, 2004). It is therefore safe to say that such a model does not fare better than other empirical wear laws and an agreement between our values for α and n and the ones of Kim et al. (1986) should not be expected.
A parameter-free mechanistic model of the adhesive wear process of rough surfaces in sliding contact Tobias Brink, Lucas Frérot, and   Video 1: Video of a smaller simulation with side length 0.25 L discretized into 512 × 512 pixels, λ l = L/128, λ u = L/16, H = 0.7, h RMS = 0.005 L, ν = 0.3,τ = 0.006 E, G = 0.25 · 10 −6 EL, where L corresponds to the side length of the surfaces used in the main text. The normal load is F N = 0.25 · 10 −3 EL 2 . The simulation was made smaller to make the features of the model easier to see. On the left a 3D view of the bottom surface is shown, the contacting areas colored according to the wear behavior: no particle detachment (yellow), a wear particle is emitted at the present moment (red), or the junction has already emitted a particle in a previous step (gray). The smaller inset at the bottom shows the same as a top view. Filename: video1-load_250e-6.mp4.
Video 2: The same as Video 1, but with a normal load of F N = 0.5 · 10 −3 EL 2 . Due to the merging of contact junctions, the wear rate is signi cantly lower than at smaller normal loads, indicating the limits of the model for high loads and small surfaces. Filename: video2-load_500e-6.mp4.