Inferences on earthquake kinematic properties from data inversion: two different approaches

Earthquake early warning (EEW) systems have been rapidly developing over the past decade. Japan Meteorological Agency (JMA) has an EEW system that was operating during the 2011 M9 Tohoku earthquake in Japan, and this increased the awareness of EEW systems around the world. While longer-time earthquake prediction still faces many challenges to be practical, the availability of shorter-time EEW opens up a new door for earthquake loss mitigation. After an earthquake fault begins rupturing, an EEW system utilizes the first few seconds of recorded seismic waveform data to quickly predict the hypocenter location, magnitude, origin time and the expected shaking intensity level around the region. This early warning information is broadcast to different sites before the strong shaking arrives. The warning lead time of such a system is short, typically a few seconds to a minute or so, and the information is uncertain. These factors limit human intervention to activate mitigation actions and this must be addressed for engineering applications of EEW. This study applies a Bayesian probabilistic approach along with machine learning techniques and decision theories from economics to improve different aspects of EEW operation, including extending it to engineering applications. Existing EEW systems are often based on a deterministic approach. Often, they assume that only a single event occurs within a short period of time, which led to many false alarms after the Tohoku earthquake in Japan. This study develops a probability-based EEW algorithm based on an existing deterministic model to extend the EEW system to the case of concurrent events, which are often observed during the aftershock sequence after a large earthquake. To overcome the challenge of uncertain information and short lead time of EEW, this study also develops an earthquake probability-based automated decision-making (ePAD) framework to make robust decision for EEW mitigation applications. A cost-benefit model that can capture the uncertainties in EEW information and the decision process is used. This approach is called the Performance-Based Earthquake Early Warning, which is based on the PEER Performance-Based Earthquake Engineering method. Use of surrogate models is suggested to improve computational efficiency. Also, new models are proposed to add the influence of lead time into the cost-benefit analysis. For example, a value of information model is used to quantify the potential value of delaying the activation of a mitigation action for a possible reduction of the uncertainty of EEW information in the next update. Two practical examples, evacuation alert and elevator control, are studied to illustrate the ePAD framework. Potential advanced EEW applications, such as the case of multiple-action decisions and the synergy of EEW and structural health monitoring systems, are also discussed.

The initial objective of Part I was to determine the nature of upper mantle discontinuities, the average velocities through the mantle, and differences between mantle structure under continents and oceans by the use of P'dP', the seismic core phase P'P' (PKPPKP) that reflects at depth d in the mantle. In order to accomplish this, it was found necessary to also investigate core phases themselves and their inferences on core structure. P'dP' at both single stations and at the LASA array in Montana indicates that the following zones are candidates for discontinuities with varying degrees of confidence: 800-950 km, weak; 630-670 km, strongest; 500-600 km, strong but interpretation in doubt; 350-415 km, fair; 280-300 km, strong, varying in depth; 100-200 km, strong, varying in depth, may be the bottom of the low-velocity zone. It is estimated that a single station cannot easily discriminate between asymmetric P'P' and P'dP' for lead times of about 30 sec from the main P'P' phase, but the LASA array reduces this uncertainty range to less than 10 sec. The problems of scatter of P'P' main-phase times, mainly due to asymmetric P'P', incorrect identification of the branch, and lack of the proper velocity structure at the velocity point, are avoided and the analysis shows that one-way travel of P waves through oceanic mantle is delayed by 0.65 to 0.95 sec relative to United States mid-continental mantle. A new P-wave velocity core model is constructed from observed times, dt/dΔ's, and relative amplitudes of P'; the observed times of SKS, SKKS, and PKiKP; and a new mantle-velocity determination by Jordan and Anderson. The new core model is smooth except for a discontinuity at the inner-core boundary determined to be at a radius of 1215 km. Short-period amplitude data do not require the inner core Q to be significantly lower than that of the outer core. Several lines of evidence show that most, if not all, of the arrivals preceding the DF branch of P' at distances shorter than 143° are due to scattering as proposed by Haddon and not due to spherically symmetric discontinuities just above the inner core as previously believed. Calculation of the travel-time distribution of scattered phases and comparison with published data show that the strongest scattering takes place at or near the core-mantle boundary close to the seismic station. In Part II, the largest events in the San Fernando earthquake series, initiated by the main shock at 14 00 41.8 GMT on February 9, 1971, were chosen for analysis from the first three months of activity, 87 events in all. The initial rupture location coincides with the lower, northernmost edge of the main north-dipping thrust fault and the aftershock distribution. The best focal mechanism fit to the main shock P-wave first motions constrains the fault plane parameters to: strike, N 67° (± 6°) W; dip, 52° (± 3°) NE; rake, 72° (67°-95°) left lateral. Focal mechanisms of the aftershocks clearly outline a downstep of the western edge of the main thrust fault surface along a northeast-trending flexure. Faulting on this downstep is left-lateral strike-slip and dominates the strain release of the aftershock series, which indicates that the downstep limited the main event rupture on the west. The main thrust fault surface dips at about 35° to the northeast at shallow depths and probably steepens to 50° below a depth of 8 km. This steep dip at depth is a characteristic of other thrust faults in the Transverse Ranges and indicates the presence at depth of laterally-varying vertical forces that are probably due to buckling or overriding that causes some upward redirection of a dominant north-south horizontal compression. Two sets of events exhibit normal dip-slip motion with shallow hypocenters and correlate with areas of ground subsidence deduced from gravity data. Several lines of evidence indicate that a horizontal compressional stress in a north or north-northwest direction was added to the stresses in the aftershock area 12 days after the main shock. After this change, events were contained in bursts along the downstep and sequencing within the bursts provides evidence for an earthquake-triggering phenomenon that propagates with speeds of 5 to 15 km/day. Seismicity before the San Fernando series and the mapped structure of the area suggest that the downstep of the main fault surface is not a localized discontinuity but is part of a zone of weakness extending from Point Dume, near Malibu, to Palmdale on the San Andreas fault. This zone is interpreted as a decoupling boundary between crustal blocks that permits them to deform separately in the prevalent crustal-shortening mode of the Transverse Ranges region.

We propose a method to obtain source parameters, most importantly the rupture velocity, of small earthquakes on the basis of stopping phases. It is assumed that a rupture begins at one focus of an ellipse and spreads radially at a constant rupture velocity on the fault surface, eventually terminating on the elliptical boundary. This model allows us to investigate a wide variety of rupture configurations, ranging from a circular to a unilaterally propagating rupture. In this model we find that two high‐frequency stopping phases, Hilbert transformations of each other, are radiated and that the difference in arrival times between the two stopping phases is dependent on the average value of the rupture velocity, the radius of the major semiaxes of the ellipse, and its ellipticity. These parameters can be estimated by a nonlinear least squares inversion method. We perform simulation tests, which allow us to detect a Hilbert transform pair whose arrival times coincide with ones predicted by an assumed source model from noisy data. These simulation tests show that our inversion method is useful. Since our method does not require an a priori assumption of the average rupture velocity, the orientation of the fault plane, and the circular crack model, it provides us with a tool for estimating accurate source dimensions as well as rupture processes of small earthquakes.

Source parameters estimated in the frequency domain for 100 selected seismic events from the Rudna copper mine, with moment magnitude ranging from 1.4 to 3.6, were collected to study their scaling relations and to compare them with the parameters estimated in the time domain. The apparent stress and static stress drop, corrected for the limited bandwidth recording, increase slightly in a similar manner with increasing seismic moment. The ratio of apparent stress to static stress drop, a measure of radiation efficiency, is practically constant and its mean value is close to 0.1. For 37 seismic events, with moment magnitude between 1.9 and 3.4, source parameters were estimated in the time domain from relative source time functions, that displayed unilateral rupture propagation, and their rupture velocity could be estimated. It ranges from 0.23 to 0.80 of shear wave velocity and is almost independent of seismic moment. The fault length, estimated from the average source pulse width and rupture velocity, is clearly dependent on seismic moment and is smaller than the source radius estimated from the corner frequency on the average by about 25 percent. There is no correlation between the values of static stress drop estimated in the frequency and time domains, but the time domain stress drop is in general similar to that estimated in the frequency domain. The apparent stress increases with increasing rupture velocity, and the ratio of apparent stress to static stress drop seems also to depend on rupture velocity.

Aftershock sequence of 19th May, 2009 Lunayyir earthquake (Mw 5.7) has been recorded by deploying seismic stations immediately after the occurrence of the main shock. Analysis of this sequence clarified that; the major part of cumulative seismic moment released during the first hours after the occurrence of the main shock; their orientation is north west (NW) parallel to the Red Sea axial trend along Najd faulting system; it is clustered at two depths; from 5-10 and 15-25 km beneath Lunayyir area; illustrated successive periods of maxima and minima; decayed rapidly following the relation of n(t)=37.28 t-0.6. Fault plane solutions for five large events indicate normal faulting is the major mechanism. These findings are in agreement with the north east (NE) (transform fault) trend runs crossing the Red Sea into Shield area where the activity has been initiated and travelled through the Lunayyir area. While the NW Najd faulting system runs parallel to the Red Sea axial trend and intersected with NE fault trend underlying Lunayyir area. Accordingly, it can be concluded that, the recent earthquake activities at Lunayyir lava flow area are related to the present-day Red Sea tectonics. Monitoring earthquake activities through installing permanent seismic network around the area is highly recommended. Key words: aftershock sequence, fault plane solution, stress, and Lunayyir.

The space-time distribution of coseismic slip of the 2011 February 21, Mw 6.2, Christchurch earthquake, New Zealand, is explored, differently from all previous studies, through a joint inversion of geodetic and strong-motion data. The geodetic data consist of both global position system (GPS), from campaign and continuous stations, and synthetic aperture radar (SAR) interferograms from two ascending satellite tracks. The strong motion data consist of 10 stations located in the Canterbury plains, these stations offering a good azimuthal coverage of the event. The kinematic rupture model for the analysed event was obtained using the parametrization and non-linear inversion scheme proposed by Delouis et al. In particular, for any subfault we explore for the local source time function (local slip history), slip direction and rupture onset time. The geometry of the fault plane used for the kinematic inversion is inferred from the analysis of the geodetic data. To validate our results we perform a resolution study for both the single and complete data sets, and an errors analysis of our final kinematic rupture model. Considering the complexity highlighted by superficial deformation data, we adopted a fault model consisting of two partially overlapping segments, with dimensions 15 × 11 and 7 × 7 km2, corresponding to different faulting types. This two-fault model, instead of single-fault model, is needed to reconstruct the complex shape of the superficial deformation data. The total seismic moment resulting from the joint inversion is 3.0 × 1025 dyne * cm (Mw = 6.2) with an average rupture velocity of 2.0 km s-1.

I determine 505 fault plane solutions from the first motions of P-waves for the background seismicity (3.0≤ M≤6.0, 1981-1991) and collect mechanisms of major earthquakes (M≥6.0, 1927-1994) from the literature in the southern California region. Then I study the seismic strain and tectonic stress fields in individual domains (ten in total) by analyzing these mechanism data. The seismic strain tensors are obtained by tensorial summation of individual seismic moment tensors. The tectonic stress tensors are determined by performing numerical inversions of the slip vector data, using Angelier's (1990) method. The findings are summarized as follows: (1) Of the 505 fault plane solutions for the 1981-1991 background seismicity, 54% are strike-slip (SF), 21 % reverse (RF), 17% normal (NF), and 8% oblique-slip faulting (OS) events. The catalog of the major earthquakes for the period 1927- 1994 also displays similar proportions of the faulting mechanisms; (2) The similarity of the focal mechanisms can be measured by a parameter, seismic consistency (Sc) introduced by Apperson (1991). It is defined as the ratio of the scalar moment of the total moment tensor to the sum of the scalar moments of individual moment tensors. In southern California, the Brawley fault (BYF) domain shows the highest Sc (0.70), whereas the White Wolf fault (WWF) domain displays the lowest Sc (0.44). Sc values in other domains vary between the above two values; (3) The depths of possible low-angle faults inferred from the fault plane solutions vary from 20 km in the Transverse ranges where N-S convergence dominates, to only 1 km in the southern Sierra Nevada fault (SSNF) domain where E-W divergence dominates. Our current data do not show the existence of a sigle unified seismically-active master detachment in the seismogenic zone; (4) The axes of the maximum principal stress, δ1. are oriented N6°E ±11°, whereas those of the maximum principal strain, e1 are oriented N5°E ±21°; (5) The strain and stress tensors are similar to each other in the Mojave (MVE), San Jacinto (SJF), Elsinore (ESF), BYF, western and eastern Transverse Ranges (WTR, ETR) domains, but dissimilar in the central Transverse Ranges (CTR), Newport-Inglewood fault (NIF), WWF, and SSNF domains. Areas with small values of Φ = (δ2 - δ3)/(δ1 - δ3) (<0.35) such as the WTR, CTR, and NIF domains are associated with more than 40% of RF events. Areas with Φ values around 0.5 such as the SJF, ETR, WWF, ESF, BYF, and MVE domains are associated with more than 47% of SF events. The SSNF domain has a large Φ (>0.65) and shows 49% of NF events, Variation of the state of stress appears to be in the Transverse Ranges where hypocenters are generally deep. Other areas show a relatively stable state of stress throughout the seismogenic depth; (6) Seismic fraction of deformation, η, is a measure of the deformation mode. It is defined as the ratio of seismic strain rate to the total deformation rate. Because of the limited seismic data, we can usually estimate the apparent instead of the real seismic fraction of deformation. Therefore, caution must be exercised in applying the values of η to evaluations of seismic potential, In southern California, there are some indications that areas in which seismic deformation nearly accounts for the total deformation are typically associated with cold and rigid batholithic rocks or high seismic velocity anomalies such as in the SJF, south central MVE, WWF, and possibly the ETR domains. However, areas with low seismic velocity anomalies are not free of earthquakes as seen, for example, in the BYF domain, which shows η = 0.6-1.0. Other domains show η < 0.4. The problem of whether the missing deformation is being released aseismically or has accumulated elastically remains to be resolved.

The paper presents a Paticles Swarm Optimization (PSO) and Simplex integrated approach to determine source parameters using InSAR data. In earthquake source parameter inversion, thereare differences in the sensitivity of different parameters to the objective function and a trade-off between the fault area and slip. We proposed that: (1) a local PSO algorithm is employed to perform a global search,(2) a 1-D Probability Desity Function (PDF) is utilized to describe the local minimum in the full parameters space, which is automatically recognized as the initial model in the Simplex algorithm inversion. The PSO/Simplex integrated approach has been successfully applied to three different earthquakes: the 2003 Bam (Iran) earthquake (strike-slip), the 1998 Zhangbei (China) earthquake (dip-slip), and a simulated event with a normal fault with right-lateral slips. All the case studies sugguest the robustness and reliability of the integrated approach.

Crosshole seismic tomography often is applied to image the velocity structure of an interwell medium. If the rocks are anisotropic, the tomographic technique must be adapted to the complex situation; otherwise, it leads to a false interpretation. We propose a nonlinear kinematic inversion method for crosshole seismic tomography in composite transversely isotropic media with known dipping symmetry axes. This method is based on a new version of the first-order traveltime perturbation equation. It directly uses the derivative of the phase velocity rather than the eigenvectors of the body-wave modes to overcome the singularity problem for application to the two quasi-shear waves. We applied an iterative nonlinear solver incorporating our kinematic ray-tracing scheme and directly compute the Jacobian matrix in an arbitrary reference medium. This reconstructs the five elastic moduli or Thomsen parameters from the first-arrival traveltimes of the three seismic body waves (qP, qSV, qSH) in strongly and weakly anisotropic media. We conducted three synthetic experiments that involve determining anisotropic parameters for a homogeneous rock, reconstructing a fault embedded in a strongly anisotropic background, and imaging a complicated four-layer model containing a small channel and a buried dipping interface. We compared results of our nonlinear inversion method with isotropic tomography and the traditional linear anisotropic inversion scheme, which showed the capability and superiority of the new scheme for crosshole tomographic imaging.

In this study, we investigate the rupture history of the 2009 April 6 (Mw 6.1) L'Aquila normal faulting earthquake by using a non-linear inversion of strong motion, GPS and DInSAR data. Both the separate and joint inversion solutions reveal a complex rupture process and a heterogeneous slip distribution. Slip is concentrated in two main asperities: a smaller shallow patch of slip located updip from the hypocentre and a second deeper and larger asperity located southeastwards along-strike direction. The key feature of the source process emerging from our inverted models concerns the rupture history, which is characterized by two distinct stages. The first stage begins with rupture nucleation and with updip propagation at relatively high (∽4.0 km s−1), but still subshear, rupture velocity. The second stage starts nearly 2.0–2.5 s after nucleation and it is characterized by the along-strike rupture propagation. The largest and deeper asperity fails during this stage of the rupture process. The rupture velocity is larger in the updip than in the along-strike direction. The updip and along-strike rupture propagation are separated in time and associated with a Mode II and a Mode III crack, respectively. The comparison between the source models inferred in this study with the Poisson ratio anomalies in the crustal volume containing the fault plane allows the interpretation of the delay in along-strike rupture propagation in terms of a structural control of the rupture history. Our results show that the L'Aquila earthquake featured a very complex rupture, with strong spatial and temporal heterogeneities suggesting a strong frictional and/or structural control of the rupture process.

As part of an effort to develop a systematic methodology for earthquake forecasting, we use a simple model of seismicity based on interacting events which may trigger a cascade of earthquakes, known as the Epidemic-Type Aftershock Sequence model (ETAS). The ETAS model is constructed on a bare (unrenormalized) Omori law, the Gutenberg-Richter law and the idea that large events trigger more numerous aftershocks. For simplicity, we do not use the information on the spatial location of earthquakes and work only in the time domain. We demonstrate the essential role played by the cascade of triggered seismicity in controlling the rate of aftershock decay as well as the overall level of seismicity in the presence of a constant external seismicity source. We offer an analytical approach to account for the yet unobserved triggered seismicity adapted to the problem of forecasting future seismic rates at varying horizons from the present. Tests presented on synthetic catalogs validate strongly the importance of taking into account all the cascades of still unobserved triggered events in order to predict correctly the future level of seismicity beyond a few minutes. We find a strong predictability if one accepts to predict only a small fraction of the large-magnitude targets. Specifically, we find a prediction gain (defined as the ratio of the fraction of predicted events over the fraction of time in alarms) equal to 21 for a fraction of alarm of 1%, a target magnitude M ≥ 6, an update time of 0.5 days between two predictions, and for realistic parameters of the ETAS model. However, the probability gains degrade fast when one attempts to predict a larger fraction of the targets. This is because a significant fraction of events remain uncorrelated from past seismicity. This delineates the fundamental limits underlying forecasting skills, stemming from an intrinsic stochastic component in these interacting triggered seismicity models. Quantitatively, the fundamental limits of predictability found here are only lower bounds of the true values corresponding to the full information on the spatial location of earthquakes.

Monitoring microseismicity is important for illuminating active faults and for improving our understanding earthquake physics. These tasks are difficult in urban areas where the SNR is poor, and the level of background seismicity is low. One example is the Newport-Inglewood fault (NIFZ), an active fault that transverses the city of Long-Beach (LB). The catalog magnitude of completeness within this area is M=2, about one order of magnitude larger than along other, less instrumented faults in southern California. Since earthquakes obey a power-law distribution according to which for each unit drop in magnitude the number of events increases by a tenfold, reducing the magnitude of completeness along the NIFZ will significantly decrease the time needed for effective monitoring. The LB and Rosecrans experiments provides a unique opportunity for studying seismicity along the NIFZ. These two array contain thousands of vertical geophones deployed for several-months periods along the NIFZ for exploration purposes. The array recordings are dominated by noise sources such as the local airport, highways, and pumping in the nearby oil fields. We utilize array processing techniques to enhance the SNR.We downward continue the recorded wave field to a depth of a few kilometers, which allows us to detect signals whose amplitude is a few percent of the average surface noise. The migrated wave field is back-projected onto a volume beneath the arrays to search for seismic events. The new catalog illuminates the fault structure beneath LB, and allows us to study the depth-dependent transition in earthquake scaling properties. Deep aseismic transients carry valuable information on the physical conditions that prevail at the roots of seismic faults. However, due the limited sensitivity of geodetic networks, details of the spatiotemporal evolution of such transients are not well resolved. To address this problem, we have developed a new technique to jointly infer the distribution of aseismic slip from seismicity and strain data. Our approach relies on Dieterich (1994)'s aftershock model to map observed changes in seismicity rates into stress changes. We apply this technique to study a three month long transient slip event on the Anza segment of the San Jacinto Fault (SJF), triggered by the remote M w 7.2, 2010 El Mayor-Cucapah (EMC) mainshock. The EMC sequence in Anza initiated with ten days of rapid (≈100 times the longterm slip rate), deep (12-17 km) slip, which migrated along the SJF strike. During the following 80 days afterslip remained stationary, thus significantly stressing a segment hosting the impending M w 5.4 Collins Valley mainshock. Remarkably, the cumulative moment due to afterslip induced by the later mainshock is about 10 times larger than the moment corresponding to the mainshock and its aftershocks. Similar to sequences of large earthquakes rupturing fault gaps, afterslip generated by the two mainshocks is spatially complementary. One interpretation is that the stress field due to afterslip early in the sequence determined the spatial extent of the late slip episode. Alternatively, the spatial distribution is the result of strong heterogeneity of frictional properties within the transition zone. Our preferred model suggests that Anza seismicity is primarily induced due to stress transfer from an aseismically slipping principal fault to adjacent subsidiary faults, and that the importance of earthquake interactions for generating seismicity is negligible.

Slow-moving landslides in high-mountain regions pose a significant natural hazard and are capable of delivering large sediment volumes to the fluvial system. Time series analysis of Interferometric Synthetic Aperture Radar (InSAR) allows us to identify unstable and potentially dangerous areas prone to landsliding, but this technique also helps quantify seasonal dynamics for predicting landslide behavior.Our study in the Eastern Cordillera of the Argentine Andes focuses on enhancing InSAR's reliability for landslide mapping. This region is characterized by moisture changes along the topographic gradient across the orogen and seasonal variability associated with the South American Summer Monsoon. We extract InSAR time series data from Sentinel-1A/B's C-band (2014-2022) and ALOS1 PALSAR's L-band (2006-2011). Tropospheric delay is caused by atmospheric turbulence and vertical stratification changes. These delays can introduce significant errors in deformation measurements, thus impacting the quality of maps portraying landslide deformation rates. To address this problem, we apply various correction techniques, ranging from spatial and temporal filtering to water-vapor estimation from an atmospheric model. Fading signal noise, another challenge caused by multi-looking and short temporal baselines in the Small Baseline Subset (SBAS) technique, additionally compromises InSAR time series accuracy. We investigate the pattern and magnitude of fading signals in landslide areas using Small Baseline Subset (SBAS) with different neighboring connections and non-linear phase inversion methods, such as the Eigenvalue Decomposition-based Maximum Likelihood (EMI), Eigenvalue Decomposition (EVD), and the Phase Triangulation Algorithm (PTA).Our research evaluates both statistical methods and Global Atmospheric Models for correcting tropospheric delays and fading signal noise. We explore statistical methods, such as double-difference filtering and corrections based on phase elevation, for different spatial windows, including individual catchments, moving windows, and adaptive window sizes. The efficiency of these methods varies with the environmental and topographic conditions in the orogen. Both stratified and turbulent components of the troposphere, along with fading signal noise, can significantly influence tropospheric delay and time series quality. In the context of the factors that influence deformation signals and the combined array of methods to obtain robust measurements, we can identify the spatial and temporal characteristics of slow-moving landslides and assess the different impacts on rate changes.

We have developed a multiple level-set method for simultaneous inversion of gravity and seismic traveltime data. The method recovers the boundaries between regions with distinct physical properties assumed constant and known, creating structurally consistent models of two subsurface properties: P-wave velocity and density. In single level-set methods, only two rock units can be considered: background and inclusion. Such methods have been applied to examples representing various geophysical scenarios, including in the context of joint inversion. In multiple level-set methods, several units can be considered, which make them far more applicable to real earth scenarios. Recently, a multiple level-set method has been proposed for inversion of magnetic data. We extend the multiple level-set formulation to joint inversion of gravity and traveltime data, improving upon previous work, and we investigate applicability of such an inversion method in ore delineation. In mineral exploration environments, traditional seismic imaging and inversion methods are challenging because of the small target size and the specific physical property contrasts involved. First-arrival seismic traveltime and gravity data complement each other, and we found that joint multiple level-set inversion is more beneficial than separate inversions, especially with limited data and slow targets. Our method is more robust than the joint inversion method based on clustering of physical properties in recovery of piecewise homogeneous models not well-constrained by the data. To justify the known property assumption, we found the degree of robustness of the multiple level-set joint inversion to uncertainties arising from incomplete knowledge of small-scale subsurface physical property variations and target composition.

Sine–Gordon modulation solutions: Application to macroscopic non-lubricant friction