How ploughing and frictional melting regulate ice-skating friction

The slipperiness of ice is well known while, for ice skating, its mechanism still needs further investigation, where the complex interactions including the thermal conduction of the skate—meltwater—ice system, the ploughing and the frictional melting of ice to the friction force are still unclear. This study presents a theoretical framework and a simplified analytical solution to unveil the friction mechanism when a curved skate sliding on ice. The theory is validated by experiments and the effects of these various factors, including the sliding velocity, the ice temperature, the supporting weight, and the geometry of the skate blade to the friction are revealed in detail. This study finds that the contribution of friction force from the ploughing deformation through skate indentation and that from the fluid friction through the shear motion of the meltwater layer is comparable with each other, which thus clarifies how the ploughing deformation of the ice substrate together with its frictional melting regulates the friction during skating.


Introduction
The slipperiness of ice is key to skaters when gliding on ice, while it causes much trouble to transportation in cold regions. So, much effort has been payed to resolve the mechanism of the slipperiness of ice for more than a century [1]. Different mechanisms have been revealed for the slipperiness of ice, such as ice surface diffusion [2][3][4], meltwater lubrication [5][6][7][8][9], mechanical ploughing [2,[10][11][12][13][14], and others [15][16][17][18][19][20], all of which are related to the specific sliding conditions. For example, when the temperature of ice is lower than about -20 °C, the slipperiness of ice will be dominated by surface diffusion, which brings an Arrheniustype equation to the coefficient of friction [2,4]. For temperature higher than -20 °C, the meltwater lubrication is widely accepted as the mechanism that defines the ultralow friction on ice, where the frictional melting [1,6,9,21], together with the pressure melting [22], and the premelting of ice [23][24][25] has been recognized as the possible explanation to the thin meltwater generation. During skating, the frictional heating is the prevailing accepted factor that creates the meltwater lubrication layer [1,8,21,[26][27][28][29]. However, due to its low viscoity, a thin water layer alone can be hardly deemed as a good lubricant [30,31]. Our recent study [32] reveals that, for a flat slider sliding on ice, a dynamic balance between the ice melting and the water squeezing regulates the thickness of the self-lubricating interstitial meltwater layer, which guarantees the low friction on ice. An analytical expression [32] is given to the coefficient of friction, which nicely captures the dependence of the friction force on other related physical variables, such as the slider velocity and its wideness, the sliding velocity and the contact pressure.
For non-flat contact on ice, which is true for ice-skating where the skate blade is always rockered   [15,21,28], both experiments and theoretical investigations have been carried to reveal the complex interplay between the friction force and the slider shape [27,29,33], the contact pressure [2,27,34], the temperature [2], and other parameters [35]. For these theoretical investigation, numerical procedure is always needed to determine the friction force on ice, while seldom of them provide a clear and selfconsistent description. For example, the well-known Frictional Algorithm using Skate Thermohydrodynamics series models (F.A.S.T.) [27,29,36,37], and their variations [9,33], deal the vertical force balance by assuming that the meltwater layer possesses a constant parabolic pressure distribution, with which its average value is equal to the hardness of ice, cannot be justified in general. Furthermore, as pointed out later by Le Berre and Pomeau [26], the condition that the growth rate of the meltwater layer thickness equals to the ice melting rate is not applicable. Although Le Berre and Pomeau [26] updated this condition by a Stefan moving boundary, their model did not provide a self-consistent description to the squeezing Poiseuille flow of the meltwater layer, which brings the thickness of the meltwater layer one order larger than others [6]. van Leeuwen [21] thoroughly investigated the friction force both in the ploughing and in the melting regimes by introducing a Bingham type constitution to ice when describing the indention dynamics of the slider. However, numerical procedure is needed in his model [21] to determine the friction force, which makes their role of those physical variables, such as the slider geometrical parameters, the supporting weight, the sliding velocity, in regulating the friction force is less clear. By treating the curved contact between the skate blade and the ice surface as a flat contact, van Dongen and Smeulders [28] built an analytical solution to determine the onset of lubrication by frictional heating and found that although the boundary friction plays a role at the tip of the skate blade, its contribution to the total friction is negligible. Based on the thermodynamic type of ice friction model developed by Makkonen [38], Spagni et al. [14] proposed an analytical model that considers the solid contact between two sliding asperities to interpret the rough contact during sliding on ice. However, a fitting procedure is needed to determine the real contact area. Based on Hertzian contact model, Böttcher et al. [5] proposed a procedure to determine www.Springer.com/journal/40544 | Friction the real contact area when determining the friction force on ice. Following the thermodynamic type of ice friction model [38], Yun et al. [6] construct a hydrodynamic model to determine the friction of a hemispherical ice on a flat solid surface, while only scaling relations has been obtained between the friction force and the contact parameters. Liefferink et al. [2] investigated the factors of ice temperature, contact pressure, and sliding speed to the slipperiness of ice and built a ploughing type friction model when the temperature of ice approaches to its melt point, where the effect of the meltwater layer lubrication is fully ignored. Based on experimental investigation, Lever et al. [18] proposed that the presence of the lubricating ice-rich slurries at discontinuous high-pressure zones other than the full contact water films determines the slippery of ice during skating. All these ice friction investigations indicate that the overall friction of skate on ice is much likely regulated by the complex interactions among these physical parameters, including the sliding speed, the ice temperature, the slider geometrical shape, the meltwater layer, the contact pressure, and the ploughing deformation. However, there is still lack a clear and self-consistent physical picture of how those parameters collectively determine the friction force on ice when the contact surface is none flat. In this study, for the curved skate blade, we aim to build a self-consistent physical model to reveal the friction mechanism when sliding on ice. We build a framework to model the friction on ice by including the ploughing deformation and the meltwater evolution. In the framework, the heat conduction in this skate-meltwater-ice system is investigated in detail and the melting rate of the ice substrate is determined accordingly. To determine the thickness evolution of the meltwater layer, the squeezing rate of the skate blade and the ploughing deformation of the ice substrate are introduced. Based on the framework, the thickness evolution of the meltwater layer and the friction force can be obtained through numerical solution. From the typical operation condition related to ice speed skating, we observe from the numerical solution that the ploughing region contributes the majority force to the sliding friction. Based on this result, an analytical solution is derived to the coefficient of friction, which is verified by the numerical solutions and validated by experiments. Based on the numerical solutions and the analytical results, the effects of the sliding velocity, the ice substrate temperature, the radius and the wideness of the skate blade, and the supporting weight to the coefficient of friction are revealed in detail. The analytical formulation clearly shows the dependence of the coefficient of friction as a function of the physical parameters related to this rounded contact sliding problem, which helps to the regulation of friction on ice on purpose.

Theories and results
In this study, the contact surface of the skate on ice is simplified as part of a cylinder, as shown in Fig. 1. Except the rockered shape at the macroscale, the surface of the skate blade and the ice substrate is treated as ideal flat surfaces. We restrict our study to the speed skating situation, where the friction speed is well above 1 m/s and the ice temperature is close to its melt point (see Table 1 for the operate condition). As has been well documented by previous researches [9,21,29,33,37], the elastic deformation of the ice substrate when supporting skaters is negligible smaller than its plastic deformation. So, the elastic deformation of the contact between the curved skate blade and the ice substrate is fully ignored. The contact state between the skate blade and the ice substrate can be divided into two parts: that is the leading region of the boundary lubricating region (not shown in Fig. 1) and the left hydrodynamical lubricating region. As has been demonstrated by van Dongen and Smeulders [28] and our previous study [32], the length of the boundary lubricating region during speed skating lies in the micrometer scale and  its contribution to the total friction force can be safely ignored. So, we will focus on the friction in the hydrodynamical lubricating region, where a thin meltwater layer exists to regulate the friction motion of the skate blade on ice.
In the hydrodynamical lubricating region, the shearing motion of the meltwater layer contributes shear resistance to the skate. The squeezing of the meltwater layer contributes radial oriented pressure to the skate, which not only supports the weight of the skater but also resists the motion of the skate. So, in order to determine the friction force, we need to determine the thickness of the meltwater layer together with its pressure field. There are three factors that regulate the thickness of the meltwater layer: the melting rate of ice, the squeezing rate of the skate, and the receding speed of the ice substrate. The melting rate of ice will be deduced from the thermal conduction analysis of the skate-meltwater-ice system. The squeezing rate of the skate is determined by the geometry of the skate and its sliding velocity. The receding speed of the ice substrate is regulated by the pressure of the meltwater layer, the strength of ice and its plastic contact deformation mechanism. Based on the thickness evolution of the meltwater layer, the classic Reynolds equation will be used to determine the pressure field of the meltwater layer. With the knowledge of the pressure field of the meltwater layer, together with the force equilibrium condition along the weight direction, the total contact length and the thickness of the meltwater layer will be revealed, which finalizes the friction force.
In Section 2.1, a general framework will be built to determine the friction force, which requires numerical procedures to determine the meltwater layer thickness and the friction force. Then, for the typical operate condition related to ice skating, the numerical solution will be presented in Section 2.2. With the observation from the numerical solution that the ploughing region contributes the majority of the friction force to the skate, an analytical theory will be proposed to the coefficient of friction on ice in Section 2.3. Finally, verification and validation of the analytical theory is done in Section 2.4 by comparison of the analytical results with the numerical solutions and the experimental results, respectively.

General framework
The thickness of the meltwater layer is regulated by the squeezing of the skate blade and the melting of the ice substrate. To determine the melting rate of ice, we need to make the heat transfer of this skate-meltwater-ice system clear. As can be checked in the following sections, the Reynolds number of the flow is no larger than 10. This is much smaller than the critical Reynolds number of ~2,100, above which the flow changes from laminar to turbulent [39]. So, during sliding, the flow of the meltwater layer belongs to laminar flow, where only the thermal diffusion regulated heat transfer will be considered. Considering the typical contact time c t between ice and skate, the typical length of the heat diffused can be given as where D is the heat diffusion coefficient and c / t L v  is the total contact time, L is the contact length, and v is the sliding speed (see Fig. 1). As the contact length is about ~0.05 m L in general and the sliding speed is well above 10 m/s for speed skating, the typical contact time c t is about  Based on the typical value of the heat penetrating length, we can simplify the heat transfer of this skatemeltwater-ice system. In general, the temperature field in the meltwater layer and in the ice substrate should be a function of position only if the system enters into the quasi-steady state; that is ( , , ) T T x y z  . As the contact length L is orders longer than the wideness of the skate blade and the heat penetrating length both for ice and meltwater, the heat transfer of this skate-meltwater-ice system can be simplified to a two-dimensional problem, where the change of these physical parameters along the sliding direction can be determined by the solution from the sectional evolution with time. Based on translation symmetry, the moving condition that correlate the contact time t and the contact length x can be given as Here the contact length x is defined in the region of [0, ] L .

Temperature field of the ice substrate
For the ice substrate, the heat penetrating length along the ice substrate is one order shorter than the wideness of the blade w (~1 mm). So, the heat transferred along the ice substrate can be simplified as the well-defined one-dimensional problem, where an initially isothermal half-space experiencing an instantaneous temperature rises at the boundary. Considering that the interface boundary of the meltwater and the ice substrate will keep moving as a result of ice melting, this resembles the well-known Stefan problem [40,41], where the temperature filed of the ice substrate can be given analytically as where s T is the initial ice temperature or the ice temperature in the far field, as shown in Fig. 2, m T is the melt point of ice, l h is the position of the meltwater-ice interface, which is related to the meltwater layer thickness and the indentation as where R is the radius of the skate blade. Here the condition of R L  has been used. From Eq. (3), the heat flux s J at the interface due to conduction into the ice substrate can be given by where i k is the heat conductivity coefficient of ice. Equation (5) indicates that the conducted heat through the ice substrate is regulated by both the contact time and the interface position. As the interface position is also a function of contact time, the heat flux s J is solely determined by the contact time.

Temperature field of the meltwater layer
To determine the temperature field of the meltwater layer, we can use the one-dimensional heat transfer configuration to describe the heat transfer problem. This can be justified as the thickness of the meltwater layer is about three orders smaller than the wideness of the blade. Under this condition, the Fourier number is given by . As can be checked with the parameters in the following section, the Fourier number of the meltwater layer is always larger than about 500 for most of the contact regions. This indicates that the unsteadiness of the temperature field that introduced by the heat generation through shear friction and the heat consumption through ice melting has been well removed by heat diffusion. So, the quasi-steady state condition can be used to Based on Eq. (8), the controlling Eq. (6) to the temperature of the meltwater layer can be given as The boundary conditions are as follows. For the lower surface of the meltwater layer, the temperature of the meltwater is the same as the melt point of ice: For the upper surface of the meltwater layer, the temperature is set to be constant as with the skate blade, see w is the temperature of the lower surface of the skate blade, where we have ignored the interfacial thermal resistance between skate and meltwater. Under these two boundary conditions, the temperature field of the meltwater layer can be given as Based on the temperature field of the meltwater layer, the heat flux transferred to the meltwater-ice interface can be given as

Temperature of lower surface of the skate blade
To determine the temperature field in the skate blade while in the y-z plane, we need to distinguish two time scales first. The first one is the typical contact time c t , which defines as the time that the skate blade cross over a fixed point, as introduced in Eq. (1). The second one is the total sliding contact time c  , which defines as the time that a skate blade gets in contact with the ice substrate. For speed skating, skaters will use either of their feet to skate on ice and the change frequency of the supporting foot determines the sliding contact time c  . In practice, , which is one order longer than the lower value of d h , the length along the z direction of the skate can be deemed as infinite. So, under quasi-steady state condition, the temperature in these far away region from the blade-meltwater interface will be constant with the temperature in the open air, which gives the far field boundary condition as T is the temperature in the skate blade,  measures the time since the skate blade contacts with the ice substrate, a T is the temperature of the blade in the far filed. In this study, a T is assumed to be constant with the open air. By ignoring the forced convective heat transfer through the free surface of the skate blade, the heat transfer in the skate blade can be deemed as a one-dimensional heat transfer problem, where the temperature change rate of the lower surface of the skate blade can be given as From Eqs. (13) and (14), the lower surface temperature of the skate blade can be given by To solve Eq. (15), the dependence of the thickness of the meltwater layer h on the sliding contact time  will be needed. As the evolution of the thickness of the meltwater layer h will be derived in detail in the following section, here we can treat it as a function of the sliding contact time  in general. Equation (15) can be solved by variables separation, which gives the solution as where the initial condition of   bs a 0 T T  has been used. So, the contribution of the initial temperature of the blade will decay exponentially, with the typical decay time 0  can be given as where b k is the thermal conductivity of the skate blade. It can be easily checked that the typical decay time 0  is orders shorter than the total sliding time c  , and the typical contact time c t . So, the surface temperature of the skate blade will always approach to a fixed temperature, which is given by Based on Eqs. (18) and (14), we demonstrate that there is no heat current to the skate blade in practice; that is b 0 J  . Although this condition has been used by those previous studies [27,29,37], this study provides a direct demonstration.

Melting rate of the ice substrate
Based on the heat flux transferred to the meltwater-ice interface (Eq. (11)) and the heat flux at the interface due to conduction into the ice substrate (Eq. (5)), the power used to melt the ice substrate can be determined, which defines the melting rate of the ice substrate as where the surface temperature of the skate blade of Eq. (18)

Thickness evolution of the meltwater layer
The thickness of the meltwater layer is controlled by the squeezing rate of the skate along the radial direction s v and the receding rate of ice r v : The squeezing of the meltwater layer is realized by indentation of the skate blade. By ignoring the elastic deformation of the skate blade, the squeezing rate of the blade can be given by where  is the angular position of the contact position, as shown in Fig. 1. As the radius of the skate blade R is orders of larger than the total contact length of L, the condition of sin   has been used.
For the receding rate of the ice substrate, we claim that the receding rate of ice is not necessarily equals to the melting rate of ice. In fact, there are three scenarios to the receding motion of the ice substrate, as shown in Fig. 3. General speaking, the receding rate of ice is controlled by the melting rate of ice 2044 Friction 11(11): 2036-2058 (2023) | https://mc03.manuscriptcentral.com/friction (given in Eq. (19)) and the plastic deformation of the ice substrate. The receding rate controlled by the plastic deformation is determined by the squeezing flow along the wideness direction of the skate blade. If there is no squeezing flow along the wideness direction, as shown in Fig. 3(a), a pure ploughing type deformation happens to the ice substrate and we call is as the ploughing region. Figure 3(b) shows the scenario when the squeezing flow along the wideness direction and the plastic deformation of the ice substrate exist together, where we call it as the mixed region. If there is no plastic deformation happens to the ice substrate, as shown in Fig. 3(c), the receding rate of the ice substrate is controlled by the melting rate of ice alone and we call it as the squeezing region.
For the ploughing region, by ignoring the squeezing flow along the wideness direction, the thickness of the meltwater layer can be determined by mass conservation as where i  is the density of ice. Taking together Eqs. (4), (19), and (22), the thickness of the meltwater layer can be determined from where the transition condition of Eq. (2) has been used and the parameter l h is given by The term i w /   in Eq. (23) originates from the density change when ice melted into water. The initial condition can be set as 0 0 x h   , where the thickness of the permelting layer of ice is ignored. The ending condition is determined by the situation that the mixed receding motion of the ice substrate happens.
For the squeezing region, the receding of ice is determined by the melting rate of ice: r m v v  . So, combining Eqs. (19)-(21) will give the control equation to the thickness of the meltwater layer The initial condition of Eq. (25) is determined by the continuity condition where the thickness and the pressure of the meltwater layer transits smoothly from the mixed region to the squeezing region. For the mixed region, the receding rate of the ice substrate is not known as a prior. In this region, as both the squeezing flow along the wideness direction and the plastic deformation of the ice substrate happens together, the thickness evolution of the meltwater will be determined by the condition that a critical pressure field of the meltwater layer has been developed, which is equals to the hardness of the ice substrate. This concept has been used by previous researches [27,29,37] and we will give it in the subsection, see Eq. (30).

Pressure of the meltwater layer
The pressure filed of the meltwater layer will be determined by the squeezing flow of the meltwater layer and the plastic deformation of the ice substrate. For the ploughing region as shown in Fig. 3(a), the pressure of the meltwater layer is equal to the hardness of ice. By ignoring the volume change during bearing the pressure load, the hardness of ice near the ice-meltwater interface will be determined by the hardness of ice in the far field with the constant www.Springer.com/journal/40544 | Friction temperature of s T . So, for the ploughing region, the average pressure of the meltwater layer over the thickness direction of the blade w ( ) p x is given by where   s H T is the hardness of ice. In this study, the hardness of ice is assumed to be a function of the ice substrate temperature only although other factors such as the loading velocity, the microstructure of ice will surely show their roles [2,42].
For the squeezing region as shown in Fig. 3(c), there is no plastic deformation happens to the ice substrate and the pressure filed of the meltwater layer will be determined by the squeezing flow via the classic Reynolds equation [43,44]. Under the condition that the contact length is much longer the wideness of the skate blade, the approximation of an infinitely short journal bearing can be used to model the squeezing flow. The simplified Reynolds equation [44,45] under this condition can be given as where ( , ) p x y is the pressure in the meltwater layer. The first term of the right hand of Eq. (27) is the result of the meltwater layer thickness variation along the contact length direction, which is determined by Eq. (25), whereas the second term is the result of squeezing rate of the skate along the radial direction, which is determined by Eq. (21). So, with the knowledge of the thickness evolution of the meltwater layer and the squeezing rate of the skate, the pressure of the meltwater layer can be determined by integrating Eq. (27) twice with the boundary condition of ( This gives the average pressure of the meltwater layer over the thickness direction of the blade w ( ) p x as Equation (29) indicates the average pressure of the meltwater layer along the sliding direction is determined by its contact position x and the thickness variation of the meltwater layer over the contact position dh/dx. For typical sliding condition, the average pressure over contact position x has been shown in Fig. 6. As can be checked from Fig. 6, the pressure given in Eq. (29) may reach to a negative value when the position x approaches to the rear point of the contact region. This negative pressure is ignored as it bounded by the cavitation pressure of the meltwater layer [46], which is much smaller than the hardness of ice.
For the mixed region, as shown in Fig. 3(b), the exist of the squeezing flow along the wideness direction implies that the classic Reynolds equation given in Eq. (27) applies, while the first term of the right hand dh/dx is undetermined. To determine this term, or the thickness evolution of the meltwater layer, the situation that the average pressure of the meltwater layer equals to the hardness of ice will be used. From Eq. (29), this gives the thickness of the meltwater layer in the mixed region as The initial condition of Eq. (30) is determined by the continuity condition where the thickness of the meltwater layer transits smoothly from the ploughing region to the mixed region; that is the ending condition of Eq. (23). The end condition of Eq. (30) is given by the initial condition of Eq. (25).

Friction force and coefficient of friction
Up to now, the thickness evolution and the average pressure of the meltwater layer is fully resolved. With the thickness evolution of the meltwater layer, the shear stress can be readily determined from . For typical sliding condition, the shear stress and the pressure of the meltwater layer has been shown in Fig. 6. Based on the pressure and the shear force in the meltwater layer, the vertical force balance can be given as where W is the weight the skate supported. This force balance condition defines the total contact length L between the skate blade and the ice substrate. Once the total contact length obtained, the friction force can be given as which gives the coefficient of friction as As there is no analytical solution to the meltwater layer thickness in general, the coefficient of friction needs to be determined numerically.

Operating condition
The analytical model proposed in Section 2.1 shows that various physical parameters will influence the skate friction force on ice. In Table 1, we summarize all these parameters and give their ranges and their typical value that we will investigate in the following sections. The hardness of the ice substrate is key to the skate blade indentation and the friction force, while its value spans in a large range in Refs. [2,42,47]. We have measured the in-situ hardness of ice in the artificial ice rink by the dropping ball method and correlated it with the ice substrate temperature (see Appendix A1 for details), as has been done by Poirier et al. [42], which gives where s T is the ice surface temperature in degrees Celsius. We remark that the measured hardness of ice disperses to a large extend even for the same ice temperature and the same ice rink, as shown in Fig. A1. This implies the hardness of ice is regulated by many other parameters, which deserves a system investigation.

Case study
Based on the typical value of these physical parameters given in Table 1, we make a case study in this section to reveal the temperature field of this sliding friction system, the thickness evolution, the pressure, and the force generation of the meltwater layer.
From literatures, the thermal conductivity of the skate blade seems plays a role in determining the friction force on ice [48][49][50] while others assumed an adiabatic boundary condition to the skate blade [27,29,37]. We check it by our analytical model and find that the adiabatic boundary condition is applicable to the skate blade. The results shown in Fig. 4(a) indicates that no matter the initial temperature of the skate blade, the surface temperature of the skate blade will approach to a stable value, which is solely determined by the melt point of ice, the thermal heat www.Springer.com/journal/40544 | Friction conductivity and the viscosity of water and the sliding velocity as given by Eq. (18). The approaching time is regulated by the initial temperature of the skate blade and its heat-conducting property. Higher temperature and larger decay time 0  as given in Eq. (17) will require a longer time for the skate blade reaches to the stable temperature. From Fig. 4(a), we can check that the time that needed for the skate blade reaches to its stable temperature lies in the microsecond and sub microsecond scale even for the SS304 blade with an initial temperature as high as 20 °C. While this time scale is at least three order shorter than the typical contact time , we can safely use the adiabatic boundary condition to the skate blade: there is no heat conduction along the skate blade.
The influence of the thermal condition of the skate blade to the meltwater layer evolution is given in Fig. 4(b). It demonstrates that higher temperature and higher thermal conductivity of the skate blade will help to the growth of the meltwater layer before the surface of the skate blade reaches to the stable temperature. However, when the surface of the skate blade reaches to the stable temperature, the thickness of the meltwater layer converges to the same value no matter the initial temperature or the heat-conducting property of the skate blade. This further demonstrates the insignificance of the thermal condition of the skate blade to the friction force generation during sliding.
Next, we determine the temperature of the ice substrate. As shown in Fig. 5 of the temperature field in the ice substrate, it clearly shows that the meltwater layer will warm the ice substrate to a certain depth, which is sliding contact length dependent. This is because warming the ice substrate is a thermal diffusion process; longer sliding length will give more time for the thermal diffusing into the ice substrate. For the typical contact length of 50 mm, the thermal penetrated length to the ice substrate is about 80 μm, which is close to the estimated value given in Eq. (1). As this thermal penetrated length is one order shorter than the wideness of the skate blade, this justifies the applicability of approximating the heat conduct in the ice substrate by a one-dimensional heat transfer process.
The thickness evolution of the meltwater layer, its averaged pressure over the wideness direction and the shear stress as a function of the contact length is determined, as shown in Fig. 6. It shows that the thickness evolution, the averaged pressure, and the shear stress of meltwater layer shows different trends in these different regions. As for the value of the meltwater layer thickness, the maximum thickness is about 1.1 μm, which is three orders smaller than the wideness of skate blade. This justifies the usage of the one-dimensional configuration to describe the heat transfer of the meltwater layer.
The ploughing region happens at the leading contact region, which characterizes a quick increase of the thickness while keeps a constant pressure. In this region, there is no squeezing flow along the wideness direction and the pressure is determined by Fig. 6 Thickness, mean pressure, and shear stress of the meltwater layer as a function of the sliding length. The mean pressure measured in MPa whereas the shear stress measured in kPa.

2048
Friction 11(11): 2036-2058 (2023) | https://mc03.manuscriptcentral.com/friction the hardness of the ice substrate, as has been given in Eq. (26). This is because, in this region, the squeezing rate of the skate blade is so fast and the thickness of the meltwater layer is so thin that the squeezing flow along the wideness direction can hardly be built before the ice substrate has been ploughed away by plastic deformation. So, all the melted ice substrate is contributed to the thickness increment of the meltwater layer.
When the thickness of the meltwater layer reaches to a certain value and the squeezing rate of the skate blade decreases to a certain value, the squeezing flow along the wideness direction can be formed, which characterizes the onset of the mixed region, where both the ploughing deformation of the ice substrate and the squeezing flow of the meltwater layer along the wideness direction happens. In this region, the thickness of the meltwater layer decreases with sliding length, which indicates the squeezing rate of the skate blade is still larger than the melting rate of the ice substrate. As can be checked, the receding rate of the ice substrate in this region is still larger than the melting rate of ice, meaning the plastic ploughing deformation still happens on the ice substrate. So, this guarantees that the average pressure of the meltwater layer has a constant value, which equals to the hardness of ice.
When approaching to the rear of the skate blade, the receding rate of the ice substrate will approach to the melting rate of the ice substrate, which indicates the onset of pure squeezing region. In this region, the classic Reynolds equation is applicable to describe the squeezing flow of the meltwater layer along the wideness direction, which gives the pressure field of the meltwater layer. From the Reynolds equation, the average pressure of the meltwater layer will decrease gradually to a negative value (not shown in Fig. 6), which is caused by a relative larger melting rate of the ice substrate compared to the squeezing rate of the skate blade. In fact, as given in Eq. (21), the squeezing rate of the skate blade is zero when x L  .
This negative pressure is insignificance as cavitation will happen to this region [46] and the Reynolds equation cannot be used to describe the flow anymore. Based on the averaged pressure of the meltwater layer along the contact length, the supporting weight and the friction force are determined and shown in Fig. 7. It shows that the ploughing region contributed a majority role to support the weight and generate the friction force, while the contribution of the squeezing region to the force generation is negligibly small (less than 3%). So, we may use the force determined from the ploughing region and the mixed region to determine the coefficient of friction. This approximated coefficient of friction is 0.00402, which compares well to the exact value of 0.00399 where the squeezing region is included. This demonstrates that, for the typical condition we considered, the contribution of the squeezing region to the coefficient of friction on ice is ignorable.

Analytical theory
In this section, we move to develop an analytical formulation to the coefficient of friction. Based on the theoretical model given in Section 2.1 and the case study given in Section 2.2, we find the contact region can be divided into three subregions and the contribution of the squeezing region to the force generation is ignorable while the majority contribution to the force generation is the ploughing region (see Fig. 7 for details). Based on this observation, we may use the solution in the ploughing region to give a suitable analytical formulation to the coefficient of friction. To do so,, we need to solve the thickness evolution of Eq. discrepancy to the thickness of the meltwater layer when compared to the numerical solution. Under this condition, the thickness of the meltwater layer can be given analytically as where the initial condition of 0 0 x h   has been used. From Eq. (35), the shear stress of the meltwater layer can be given as The average pressure will be constant with the hardness of the ice substrate.
With the knowledge of the shear stress and the pressure of the meltwater layer, the vertical force balance given in Eq. (31) can be solved, which gives the total contact length: where only the leading term related to / L R is reserved. The contribution of the shear stress to the weight support is ignored because of its minor contribution.
The friction force now turns to where only the leading term related to / L R is reserved again. Finally, the coefficient of friction can be given analytically from Eqs. (37) and (38) where m s T T T    is the temperature difference between the melt point and the ice substrate in the far field, the constants 1 C and 2 C are determined by the physical properties of water and ice, which are given by Equation (39) gives the analytical dependence of the coefficient of friction as a function of the supporting weight, the radius and the wideness of the skate blade, the sliding velocity, and the temperature difference and the hardness of ice, while with no fitting parameter. It shows that the contribution to the coefficient of friction can be divided into two parts. The first part corresponds to the normal ploughing force, which is caused by the ploughing motion of the curved skate blade on the flat ice substrate. Liefferink et al. [2] has used this ploughing model to correlate the effect of the applied loading on the coefficient of friction in their ball-on-ice configuration when the temperature of ice closes to the melt point and the sliding velocity is low. The second part corresponds to the tangential fluid friction force, which is caused by the shear resist of the meltwater layer. This part applies when a meltwater water layer can be formed and it will scale with 1

Theory verification and validation
In this section, we make a verification to the analytical formula through a system comparison with the numerical solutions during revealing the effect of these various operate conditions to the coefficient of friction. The theory will be validated at the end of this section by experiments.

Ice substrate temperature effect
The effect of the ice substrate temperature to the coefficient of friction and the total contact length is given in Fig. 9. It shows from Fig. 9(a) that the analytical solution could capture the effect of the ice substrate temperature to the friction force. The analytical theory will give a slightly larger value to the coefficient of friction when the ice substrate temperature is high and it will give a bit lower value if the ice substrate temperature is low. The maximum discrepancy between them happens to the lower boundary of the temperature we considered, which is about 8.1%. Figure 9(a) also shows that increasing the ice substrate temperature will first decrease the value of the coefficient of friction and then increase its value. That is to say there is an optimal ice substrate temperature that the friction force reaches to its minimum. The optimal ice substrate temperature is about -4~-6 °C for the parameters we concerned, which is close to the experimental measurements [51]. When the ice substrate temperature is lower than the optimal one, decreasing the ice substrate temperature will introduce less ploughing type friction while add more fluid friction, which is caused by the reduced meltwater layer thickness. When the ice substrate temperature is higher than the optimal one, increasing its value will cause much ploughing type friction. This is caused by the reduced hardness of ice and the increased contact length between skate and ice when the ice substrate temperature approaches to its melt point. As shown in Fig. 9(b), the total contact length will increase about 1.5 folds if the ice temperature approaches to the melt point. This large increment of the total contact length explains the increase of the friction force.

Sliding velocity effect
In most of the time, the sliding velocity determines the performance of the skater. In Fig. 10, we determine the coefficient of friction as a function of the sliding velocity. Both the numerical solution and the analytical results indicate that there is an optimal velocity that the coefficient of friction reaches to its minimum, where the energy efficiency of the sliding is best. For skating under typical condition as given in Table 1, the optimal velocity is about 6.6 m/s from the numerical solution, which compares to the value of 8.4 m/s that from the analytical solution. The coefficient of friction as a function of the sliding velocity matches nicely with each other between the numerical solution and the analytical results if the sliding velocity lies in the range of 8-15 m/s. Otherwise, there is a certain discrepancy between them, which is no more than 8.9% in the sliding velocity range we considered in this study.
www.Springer.com/journal/40544 | Friction Fig. 10 Coefficient of friction as a function of sliding velocity. Figure 10 also shows that below the optimal velocity, the coefficient of friction will increase quickly if the sliding velocity decreases and especially when the velocity approaches to 0. On the other hand, if the velocity is above the optimal value, the coefficient of friction will increase slowly with increasing sliding velocity. As the friction caused by ploughing is almost constant for these various sliding velocities, the change of the coefficient of friction is mainly attributed to the fluid friction. When the sliding velocity approaches to zero, the thickness of the meltwater layer will decrease quickly as there is little frictional heat generated by the friction motion. So, a larger shear stress exists in the meltwater layer, which brings the higher coefficient of friction at the low sliding velocity.

Supporting weight effect
The supporting weight of the skate blade is not only determined by the weight of the skater, but also regulated by the loading condition. In Fig. 11, we considered the weight range spans from 0.5-fold to 5-fold of the typical weight, which are corresponding to the condition of 1 g while supported by two feet equally and the condition of 5 g while supported by one foot, respectively. Under these conditions, the analytical formulation matches nicely with the numerical solution both for the coefficient of friction and the total contact length. The maximum discrepancy for the coefficient of friction happens at the lower weight condition, which is 3.1%, while that for the total contact length happens at the higher weight condition, which is 1.4%. Furthermore, under the reference operation condition, both the analytical theory and the numerical solution indicates that there is a critical weight of about 500 N that the coefficient of friction reaches to its minimum.
As shown in Fig. 11(a), increasing the supporting weight of the skate blade will increase the coefficient of friction significantly as compared with the change of the substrate temperature or the sliding velocity, as shown in Figs. 9(a) and 10. This is caused by the quick increase of the ploughing force that the rounded skate blade encountered when supporting the heavy load. In fact, as shown in Fig. 11(b), the total contact length will increase almost linearly with the supporting weight, which indicates a linear increment of the ploughing depth and the ploughing force. This quick increase of the coefficient of friction as a function of the supporting weight has also observed in the ball-on-ice configuration [2], indicating the importance of the ploughing force to the total friction on ice.

Skate blade geometrical effect
For the rounded skate blade, these two parameters of the radius and the wideness determine the contact state of the skate on ice and the friction force. Figure 12 shows the map of the coefficient of friction as a function of the radius of the skate blade and its wideness while both from the numerical solutions and from the analytical results. It clearly shows that the results determined from the analytical solutions match well with the numerical solutions in general, with the maximum discrepancy accounts to 4.5%.
From Fig. 12, a smaller coefficient of friction can be attained if the radius and the wideness of the skate blade are larger, although the wideness of the skate blade contributes little to the coefficient of friction if the radius of the skate blade is sufficient large. The coefficient of friction will increase quickly if the radius of the skate blade approaches to a smaller value of several meters. Under this condition, the skate blade will indent significantly into the ice substrate, which introduces a quick increase of the ploughing type friction to the skate blade. This explains why the measured coefficient of friction from the ball-on-ice configuration [1] is always much larger than that of the skate-on-ice. On the other hand, if the radius of the skate blade is large, the ploughing force along the sliding direction will be largely relieved. This explains why the skate blade used by the speed skating has a radius over 20 m in general.

Comparison with experiments
In order to validate our physical model by experiment, we conducted friction force measurement on the artificial ice rink by a homemade vehicle, with details given in the Appendix A2. Table 2 gives all the operate conditions conducted by the homemade vehicle, with the supporting weight ranges from 160 to 330 N, the sliding velocity ranges from 0.75 to 5.03 m/s, the ice substrate temperature ranges from -4.14 to -6.61 °C. The in-situ hardness of the ice substrate is also measured (see Appendix A1 for details) and included in Table 2. The radius and the wideness of the skate blade is 8.15 m and 1.2 mm, respectively.
According to these operate conditions, the coefficient of friction from the experimental measurement, from the numerical solutions and from the analytical results are given in Fig. 13. A comparison between them shows that the theory developed in this study can capture the experiments nicely in general. Furthermore, the analytical formulation can always provide a meaningful forecast for these various operate conditions given in Table 2. For example, the dependence of the coefficient of friction on the sliding velocity can be obtained from the operate conditions 6-9, where  A close comparison between these different operate conditions shows that under operation counts of 5 and 11, there is a remarkable deviation (up to 22%) between the experiments and the theoretical results. This may be caused by the dispersity of the hardness of the ice substrate. As has been mentioned in the former section, the hardness of ice will be influenced by many factors. Even on the same ice rink with the same temperature and the same ice-making technology, the hardness of ice would vary as much as 74% of its average value (see Fig. A1). The change of the hardness of the ice substrate will significantly influence the indentation and the coefficient of friction. A comprehensive investigation through experiments to the friction force and the comparison with theory are under conduction, which will be published in the near future.

Discussion
In general, the skate friction on ice is a three-dimensional multiscale problem that spans from molecule scale (the permelting layer of ice) to tens of meter scale (the radius of the skate blade). So, developing a physical model that can precisely include all the information in these different scales possesses much difficult. In doing so, most of these models are based on highly simplified assumptions while without proper justification. For example, the F.A.S.T. series models [27,29,36,37] describe the sequeezing of the meltwater layer by assuming a constant pressure exist between the skate blade and the ice substrate, which fully ignored the plouging region and the squeezing region as shown in Fig. 3. These models also ignored the Stefan moving boundary condition when determining the heat flux to the ice substrate. In the model proposed by van Leeuwen [21,33], although the Stefan moving boundary condition is introduced, the assumption that half of the frictional heat are diffused away through the skate cannot be justified. In fact, based on the thoroughly analysis of the thermal conduction of this blade-meltwater-ice system, we demonstrate that the heat dissipated by the skate blade happens at the very early stage after contact, which can be ignored in practice, as shown in Fig. 4. This study proposes the theoretical framework to model the friction force on ice with certain assumptions, which have been verified carefully by the numerical solutions. To make clear the effect of all these physical parameters in regulating the friction on ice, we pursue an approximate while analytical solution to the coefficient of friction on ice by observing that the ploughing region contributes the majority force to the friction and the load support, as shown in Fig. 7. The analytical formulation is verified by the numerical solution and both of which are compared with the experimental measurements, demonstrating the suitability of the formulated theory in forecasting the friction force on ice. Our theory indicates that the thermal conduction of the skate blade shows little effect in regulating the friction force as has been observed by Liefferink et al. [2], although early researches proposed that its role is observable [49,50]. For the typical operate condition related to skating, there is a critical value for the sliding velocity and the ice substrate temperature that the coefficient of friction can reach its minimum if the supporting weight and the geometry of the skate blade are fixed. Both of these results are consistent with the trends observed from the experimental results [1,48,[50][51][52], although we cannot make a direct comparison between them as a result of the missing parameters, especially the missing information for the hardness of ice. So, we measure the coefficient of friction on ice by a homemade vehicle, together with all the related parameters needed, to validate the theory. With the limited experimental data for the sliding condition, we find the match between the theory and the measurements is fairly well in general (Fig. 13), although certain discrepancy does exist. This discrepancy is assumed to be the uncertainty encountered in experiment, which mainly attributes to the nonhomogeneous of the ice hardness. As the hardness of ice is key to determine the indention of the skate on ice, its contribution to the coefficient of friction is of great importance. This calls for the importance of conducting a system investigation to the hardness of ice, which is now under investigation by our group.
Both the theory and the experimental measurement show that the coefficient of friction for skating can be as low as 0.003. Why ice is so slippery? This study indicates that the ploughing force along the sliding direction and the shear force of the meltwater layer determines the coefficient of friction on ice. For low speeding friction while under relative low ice temperature, Liefferink et al. [2] indicated that, in their ball-on-ice configuration, the ploughing force alone can be used to determine the friction force. While the ploughing type friction is consistent with our friction model, we show the ploughing alone cannot provide the meaningful coefficient of friction if a meltwater layer do exist between ice and slider. The ploughing alone treatment assumes that there is no friction between ice and skate except ploughing even relative shear motion happens between slider and ice. This is applicable only if the shear friction between ice and slider is negligibly small when compared with the ploughing force. Although this can be justified in their study [2], it is insufficient to describe the friction if a meltwater layer does exist between ice and slider, as the propose of this study. This study indicates that the ploughing alone cannot provide the meaningful coefficient of friction. From Figs. 9-11, it clearly shows that the shear force caused by the lubrication layer contributes at least as an equal part as that of the ploughing force. This does not mean that the exist of the meltwater layer increases the friction force on ice. On the contrary, it acts as a good lubricant that separated the direct contact between the skate blade and the ice substrate, which resembles the meltwater layer during the flat contact [32]. In fact, as a result of the sufficient large radius of the skate blade compared with the indentation contact size, together with the meltwater layer regulated lubrication, the skate blade can reach the sufficient low coefficient of friction.
Both from the literatures (see Table A1) and from our measurement, the coefficient of friction for skating is always lower than about 0.01. This is remarkable if we note that other shape of sliders and other measurement procedures could bring a coefficient of friction to the value up to 0.7 [1,2], which is almost two orders larger than that of the skating. To explain the robustness of this ultralow coefficient of friction related to skating, we can resort to the analytical formula given in Eq. (39). The first term of the coefficient of friction of Eq. (39) is given by the ploughing type friction. Considering the typical ice hardness of 10 MPa, the typical weight of the skater of 700 N, the wideness of the skate blade of 1 mm, and the radius of the skate of 10 m, this term amounts to 0.0035. The second term is given by the shear friction of the meltwater layer. Considering the constants 1 C and 2 C given in Eq. (40), which can be determined to be 1.1513×10 3 °C·Pa·s 3/2 /m and 0.8952 °C 2 ·s 2 /m 2 , respectively, and the sliding velocity of 10 m/s and the ice temperature of -5 °C, this term amounts to 0.0024. Token together these two terms, it shows that the typical coefficient of friction is 0.0059, which is well below 0.01. This implies the ultralow www.Springer.com/journal/40544 | Friction friction on ice during skating is regulated by the geometrical properties of the skate and the meltwater layer evolution. While the meltwater layer formation is the outcome of the frictional heat, the morphology of the skate blade attributes to the robustness of the ultralow coefficient of friction related to skating.
Although this study has systematically investigate the effect of the sliding velocity, the radius, the wideness, the thermal properties and the supporting weight of the skate blade, and the temperature and the hardness of ice to the friction on ice, there are still other key factors that are not included now, such as the wettability [16] and the surface roughness [12,14,16,18,53] of the skate blade, the relative humidity of the open air [54], and the elastic deformation of the system. Based on our theoretical model, we remark that the wettability and the surface roughness of the skate blade will show their role to the meltwater layer evolution while have little contribution to the ploughing force if the plastic deformation dominants the deformation of the ice substrate. The existence of roughness on the skate blade indicates the true contact happens at these asperities, which implies the evolution of the meltwater layer will be determined by considering all these asperities together. For rounded contact and flat contact, this study and our previous study [32] has revealed the details of the meltwater evolution. So, in order to reveal the friction on ice with contact roughness, we need to extend the study to a broader shape of contacts, such as ball-on-ice configuration. Furthermore, the contact roughness together with the wettability of the skate blade will influence the friction force through capillary bridge, which is significant if the roughness lies in micro-and nanoscale. So, the effect of the complex interaction between roughness and wettability to the friction on ice deserves further investigation.

Conclusions
In conclusion, this study presents a physical model to reveal the mechanism of the skate friction on ice with rounded contact. The thermal conduction of the skate-meltwater-ice system is investigated in detail and from which the melting rate of the ice substrate is determined. Together with this melting rate, the squeezing rate of the skate blade and the ploughing deformation of the ice substrate, a theoretical framework is formulated to reveal the thickness evolution of the meltwater layer and the friction force generation. Observed from the numerical solution that the ploughing force is found to contribute the majority force to the friction, an analytical formulation is developed to the coefficient of friction, which is verified by the numerical solutions and validated by experiments. The analytical model reveals that the ploughing force along the sliding direction and the fluid friction caused by the meltwater layer together contribute the friction force. For ice skating, their contribution is comparable with each other. Based on the theoretical model, the effects of these factors, including the sliding velocity, the ice temperature, the radius and the wideness of the skate blade, and the supporting weight, to the coefficient of friction are revealed in detail. The formulated explicit relation between the coefficient of friction and the related physical parameters will provide valuable insight to the regulation of the friction on ice. The temperature of the ice substrate is measured by thermocouple thermometer. Figure A1 gives a collection of the measure hardness of ice and the corresponding ice temperature, which is fitted by a linear regression equation.

A2 Friction measurement
To measure the coefficient of friction on ice, we designed a homemade vehicle that is drove by high power fans. Two identical skate blades are mounted on the vehicle tightly. The weight of the vehicle is regulated by carrying the dumbbells with different weight on propose. The accelerometer is placed on the vehicle to record the acceleration. The velocity of the vehicle is determined by the high-speed photography method. First, the vehicle is accelerated by the fans to a certain velocity. Then the fans are turned off and the vehicle is left slide freely on ice. The accelerometer on the vehicle records the deceleration of the vehicle. The air drag of the vehicle is determined by the wind tunnel test under the same velocity. Based on the deceleration of the vehicle on ice and the air drag, the coefficient of friction on ice can be determined.  The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

A3 A literature collection of the experimental measured coefficient of friction related to skating
To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.