Propulsion Calculated by Force and Displacement of Center of Mass in Treadmill Cross-Country Skiing

This study evaluated two approaches for estimating the total propulsive force on a skier’s center of mass (COM) with double-poling (DP) and V2-skating (V2) skiing techniques. We also assessed the accuracy and the stability of each approach by changing the speed and the incline of the treadmill. A total of 10 cross-country skiers participated in this study. Force measurement bindings, pole force sensors, and an eight-camera Vicon system were used for data collection. The coefficient of multiple correlation (CMC) was calculated to evaluate the similarity between the force curves. Mean absolute force differences between the estimated values and the reference value were computed to evaluate the accuracy of each approach. In both DP and V2 techniques, the force–time curves of the forward component of the translational force were similar to the reference value (CMC: 0.832–0.936). The similarity between the force and time curves of the forward component of the ground reaction force (GRF) and the reference value was, however, greater (CMC: 0.879–0.955). Both approaches can estimate the trend of the force–time curve of the propulsive force properly. An approach by calculating the forward component of GRF is a more appropriate method due to a better accuracy.


Introduction
Forces acting on a skier's center of mass (COM) in a forward direction are propulsive forces, which are the primary mechanical determinants of an cross-country (XC) skier's performance [1]. The position of skier's COM can be obtained by using the marker-based motion capture system with a segmental method [2]. Thus, forces acting on a skier's COM can be obtained by multiplying COM acceleration with the total mass of the skier, and this will indicate how athletes overcome resistive forces. However, the contribution of single pole and leg thrusts could not be revealed. Therefore, it is essential to compute forces acting on the COM from the ground reaction forces (GRFs) generated from skis and poles, separately.
Except for estimating the propulsive force with the forward acceleration of COM and the total mass, other approaches have been developed. One approach is to estimate the propulsive force as the forward-directed horizontal component of the three-dimensional (3D) GRFs from both skis and poles that act on a skier (F net ) [3]. The roller skis [3][4][5], skis [6,7], and poles [3,6,7] equipped with force sensors have been used to measure the forces generated from skis and poles. Combined with the pole angle, ski angle, ski-edging angle, and the incline of the track or the treadmill, the propulsive force from skis and poles can be specified [1,3,8]. Therefore, questions related to the propulsive force, including the contribution of skis and poles in different techniques [3,9,10], and the comparison of different techniques [11,12] have been addressed. Another approach, demonstrated by Göpfert et al. [6], is to estimate the propulsive force with the forward component of translational force (F pro ). The translational force was modeled as the component of the 3D resultant GRFs that acts in the direction from the point of force application (PFA) to a skier's COM [6], and calculated by projecting the GRFs to the line defined by the COM and PFA.
The propulsive forces obtained with the two mentioned approaches (F net and F pro ) have been compared to the propulsive force calculated with COM acceleration from a motion analysis system (F) in [6]. As using the segmental method has been shown to be suitable for estimating the position of the COM in sports [2], F was considered as the reference value. The results indicated that the force-time curves of F net and F pro all showed high similarity when compared to the force-time curves of F during the leg skating pushoffs on snow. F net overestimated F, and F pro was found to be a more appropriate approach to estimate F during leg skating push-offs [6]. However, whether F net and F pro could be used to estimate F, and which one is more accurate in other techniques, are still unknow. As XC skiing is a sport whose competition and training are normally performed on varying track topography and speed, whether F net and F pro could work steadily when estimating F at different terrain with different speeds need further investigation.
Therefore, the first aim of the present study was to obtain the force-time curves of F net and F pro with different skiing techniques and evaluate which can estimate the force-time curves of F better. As the use and importance of double poling (DP) and V2 skating (V2) as main techniques in XC skiing have increased for the past few years [12][13][14], DP and V2 skating techniques will be performed in this study. The second aim is to investigate which approach is more accurate when estimating F. Another aim is to explore the stability of the approaches to calculate F net and F pro by changing the speed and incline of the treadmill. We hypothesized that the force-time curves of F net and F pro all give comparable shape with F in both techniques [6]. We also hypothesized that F net would give a considerable overestimation, and F pro would be more accurate than F net , when estimating F in both DP and V2 techniques [6]. We further hypothesized that the approaches to calculate the F net and F pro would not be affected by the speed and incline of the treadmill in both techniques.

Participants
A total of 10 experienced male skiers (age: 29.4 ± 7.9 years; height: 181.4 ± 5.7 cm; weight: 77.9 ± 8.9 kg) who were familiar with treadmill roller skiing volunteered to participate in this study. The experimental protocol and all methods used in this study were approved by the Ethics Committee of the University of Jyväskylä. All participants provided written informed consent before the measurement and were free to withdraw from the experiments at any point.

Protocol
The anthropometric parameters needed for motion analysis (e.g., bilateral leg length, knee width, ankle width, shoulder offset, elbow width, and hand thickness) were measured first, and passive reflective markers were attached to the participants and equipment. Once the preparations were made, participants completed a 10-15 min warm-up roller skiing on the treadmill. Next, calibration was performed with the skier in a standing position and the treadmill at a 0 • incline. Participants then performed the DP technique at five speeds (13,15,17,19, and 21 km/h) on a 2 • incline. The comfortable pole length for the DP technique was 1.56 ± 0.06 m. After the trials with varying speeds at a 2 • incline, the DP technique was performed at three inclines (3 • , 4 • , and 5 • ) with a speed of 10 km/h. There was a 1 min rest between each speed and incline. When participants finished performing the DP technique, the pole length was adjusted to a comfortable length for the V2 skating technique (which was 1.63 ± 0.03 m in this study). The participants were given a short rest period while adjusting the pole length. The participants then performed the V2 technique on the treadmill. The protocol for speed and incline change was the same as during the DP test.

Data Collection
An eight-video-camera motion capture system (Vicon, Oxford, UK) and NEXUS 2.8.1 software were used to collect and record the 3D trajectories of reflective markers at a sampling rate of 150 Hz. The global coordinate system (GCS) was defined using the right-hand rule when the incline of the treadmill was 0 • . The X-axis of the GCS was defined as the direction from side to side across the treadmill. The Y-axis of the GCS was the longitudinal axis of the treadmill. The Z-axis of the GCS was perpendicular to the ground, pointing upward. The GCS was calibrated according to Vicon's specifications. A total of 58 passive reflective markers were used in this current study: 43 passive reflective markers were attached to the participants' bodies, and 15 markers were attached to the equipment, including both skis (3 each), both poles (3 each), and the treadmill (3). Anthropometric measurements and the placement of markers on the participants' bodies were conducted according to the XC model [6] used in previous studies. Measurements were performed on a motorized treadmill with a belt surface of 2.7 m wide and 3.5 m long (Rodby Innovation AB, Vänge, Sweden). The same pair of roller skis were used for both techniques (Marwe, SKATING 620 XC, wheel no. 0), with a resistance friction coefficient of µ = 0.025 measured before the measurement (Appendix A).
Two custom-made pole force sensors (VTT MIKES, Technical Research Centre of Finland Ltd., Kajaani, Finland, Figure 1a) were used to measure the axial GRF from the poles. Two custom-made two-dimensional (2D) force measurement bindings (Neuromuscular Research Centre, University of Jyväskylä, Jyväskylä, Finland, Figure 1b) [15] were mounted on roller skis to measure the leg forces generated from roller skis. Both pole and leg forces were collected synchronously with the Coachtech online measurement and feedback system (Neuromuscular Research Centre, University of Jyväskylä, Jyväskylä, Finland) at a sample rate of 400 Hz. The force measurement bindings measured the vertical (F skiz ) and mediolateral (F skix ) forces and were calibrated before the measurement [15]. A trigger signal was sent from the Coachtech [16] to the motion capture system to mark the start of the force capture. The nodes for the pole force sensors and force measurement bindings were used to supply power and transmit data.

Data Reduction
Marker labeling and COM calculations were performed using NEXUS 2.8.1 software. The raw 3D trajectories of all reflective markers and the acceleration of COM were lowpass filtered (fourth-order, zero-lag, and Butterworth filter) with a cutoff frequency of 11.3 Hz [17]. The XC model [6], which contained the head, thorax, abdomen and pelvis, upper arms, hands, thighs, shanks, feet, skis, and poles, was used to calculate the whole-  Figure 2. The segmental anthropometric data were taken from Dempster's study as described in Selbie et al. [18]. Force data were low-pass filtered (eighth-order, zero-lag, and Butterworth filter) with a cutoff frequency of 15 Hz [19]. Data filtering and parameter calculations were performed using MATLAB R2018a (MathWorks, Natick, MA, USA).    i . The PFA f and PFA r were the points of force application of the front and rear sensors, respectively. The distance between Ski_2 and PFA f was m, and the distance between PFA f and PFA r was n. markers on the roller ski ( Figure 3). The transformation from the roller ski system to the GCS is given by  where F x , F y , and F z are the components of forces generated from legs ( −→ F r ) in the GCS. R [18] is the rotation matrix from FCS to GCS.
The PFA is needed for calculating the translational force introduced by Göpfert et al. [6]. The displacement of the PFA along the binding (PFA ski ) was calculated from the force distribution between the front and rear sensors of the binding [20] over time. The PFA for each part of the binding (PFA f and PFA r ) was defined as the center of each sensor. The distance between the marker Ski_2 and PFA f (m), and the distance between PFA f and PFA r (n) were measured before the measurement. The displacements of PFA f and PFA r in GCS were obtained by moving the midpoint of Ski_2 and Ski_3 along the opposite direction of → j . The moving distances were m and m + n, respectively. The mediolateral sway of PFA ski on ski binding was not considered in this study. Thus, the PFA ski moved between the PFA f and PFA r ( Figure 3).
The measured axial pole forces ( −→ F p ) were considered the GRFs acting along the pole from the tip to the top of the pole and expressed that way in the GCS. The magnitude of −→ F p was collected using a pole force sensor. The direction of −→ F p was defined using the reflective markers that were attached to the pole. The PFA of poles (PFA p ) was defined as the intersection of the plane of the treadmill and the long axis of the pole. The plane of the treadmill was defined using the three markers attached to the treadmill.

The Reference Force, the Total Resultant Force, and the Translational Force
As using the segmental method has been shown to be suitable for estimating the position of the COM in sports [2], forces calculated by COM acceleration ( → a ) multiplied the total mass of the subject, and the equipment was the reference force ( → F ) in this study. One approach to estimate forces acting on skier's COM is to calculate the total resultant force ( −→ F net ) without considering the position of COM. −→ F net is calculated as where → G is the gravitational force of each participant and all the equipment. − −−− → F friction is the frictional force between the roller ski and the treadmill, which was directed along the path of the ski motion, and the magnitude was computed by multiplying µ with F skiz .
Another approach to estimate forces acting on skier's COM is to calculate the total translational force (the mechanical principle of translational force, see Appendix B). The translational force is the share of the resultant GRF acting in the direction from PFA to COM. The translational force from skis ( Figure 4) is the share of ski GRF ( −→ F r ) acting in the direction defined from PFA ski to COM and is calculated from where → u is the unit vector determined from PFA ski to COM. The translational force from  where v ⃗⃗ is the unit vector determined from PFAp to COM. The total translational force (F pro ⃗⃗⃗⃗⃗⃗⃗⃗⃗ ) is the sum of the translational force from the legs, poles, and the resistance. Thus, F pro ⃗⃗⃗⃗⃗⃗⃗⃗⃗ can be computed as As forces acting on skier's center of mass (COM) in forward direction are the propulsive forces, the Y component of F ⃗⃗ , F net ⃗⃗⃗⃗⃗⃗⃗⃗ , and F pro ⃗⃗⃗⃗⃗⃗⃗⃗⃗ (F, Fnet, and Fpro) was compared and analyzed in the present study.

Figure 4.
Diagram of force decomposition from skis. F r is the resultant force generated from legs. F tS is the translational component, which went through the COM. F pS represents the propulsion generated from legs in the forward direction.

Cycle Definition and Analyzed Parameters
A total of 10 consecutive poling phases for each DP technique trial and 10 consecutive kicking phases (5 left ski kicking and 5 right ski kicking) for each V2 technique trial were analyzed. The poling phase was defined as the period from the start of the pole ground contact to the end of the pole ground contact (Figure 5a). The kicking phase was defined as the ski force minima until the end of ground contact [7] (Figure 5b). The forces of skis and poles from both the left and right sides were included while calculating the total propulsion in both techniques.
The positive square root of the adjusted coefficient of multiple determination, which is the adjusted coefficient of multiple correlation (CMC, 0 < CMC < 1) [21,22], was calculated for evaluating the similarity of force-time curves. One comparison was between Fnet and F force-time curves. The similarity between Fnet and F was represented by CMCnet. The mean force difference and mean absolute force difference between Fnet and F were M F net −F and M |F net −F| . Another comparison was between Fpro and F force-time curves. The similarity between Fpro and F was represented by CMCpro. The mean force difference and mean absolute force difference between Fpro and F were M F pro −F and M |F pro −F| . The mean force differences and mean absolute force differences were computed over force curves averaged over 10 force-producing phases. The mean force differences, which are M F net −F and M F pro −F , were calculated to provide descriptive statistics only. The forces in this study were presented as values relative to body weight (%BW). −→ F net , and −→ F pro (F, F net , and F pro ) was compared and analyzed in the present study.

Cycle Definition and Analyzed Parameters
A total of 10 consecutive poling phases for each DP technique trial and 10 consecutive kicking phases (5 left ski kicking and 5 right ski kicking) for each V2 technique trial were analyzed. The poling phase was defined as the period from the start of the pole ground contact to the end of the pole ground contact (Figure 5a). The kicking phase was defined as the ski force minima until the end of ground contact [7] (Figure 5b). The forces of skis and poles from both the left and right sides were included while calculating the total propulsion in both techniques.

Statistical Analyses
A two-way mixed factorial ANOVA was performed. The dependent variables of the analyses were (1) CMCs and (2) mean absolute force differences. The independent varia- The positive square root of the adjusted coefficient of multiple determination, which is the adjusted coefficient of multiple correlation (CMC, 0 < CMC < 1) [21,22], was calculated for evaluating the similarity of force-time curves. One comparison was between F net and F force-time curves. The similarity between F net and F was represented by CMC net . The mean force difference and mean absolute force difference between F net and F were M F net −F and M |F net −F| . Another comparison was between F pro and F force-time curves. The similarity between F pro and F was represented by CMC pro . The mean force difference and mean absolute force difference between F pro and F were M F pro −F and M |F pro −F| . The mean force differences and mean absolute force differences were computed over force curves averaged over 10 force-producing phases. The mean force differences, which are M F net −F and M F pro −F , were calculated to provide descriptive statistics only. The forces in this study were presented as values relative to body weight (%BW).

Statistical Analyses
A two-way mixed factorial ANOVA was performed. The dependent variables of the analyses were (1) CMCs and (2) mean absolute force differences. The independent variables were the speed (or the incline) of the treadmill and the comparisons (i.e., comparison between F net and F and comparison between F pro and F). The speed (or the incline) of the treadmill was treated as the within-subject factor, and the comparison pair was treated as the between-subject factor. The EMMEANS subcommand with the Bonferroni adjustment in SPSS was used to perform the pairwise comparisons of the dependent variable when interactions were detected [23], and the effect size ( p η 2 ) was calculated for further evaluation. The level of statistical significance was set at α = 0.05. All data are presented as mean ± standard deviation (SD). Data analyses were conducted using version 23.0 of the SPSS program package for statistical analysis (SPSS Inc., Chicago, IL, USA).

Results
The force-time curves of F, F net , and F pro for the DP and the V2 techniques are shown in Figure 6. The interaction effect (comparison * speed) was significant on CMC with the DP technique (p = 0.038), but not with the V2 technique (p = 0.988). CMC pro did not differ from CMC net at any speed in the DP technique (p ≥ 0.106, Table 1). With the V2 technique, the overall CMC pro was about 5% lower than CMC net (p = 0.011, Table 1). The interaction effect (comparison * incline) was not significant on CMC in the DP technique (p = 0.620) but was significant in the V2 technique (p = 0.042). In the DP technique, the main effect of comparison on CMC was not significant (p = 0.218, Table 1). In the V2 technique, CMC net was significantly greater than CMC pro at 3 • (p = 0.042, Table 1).

21, 21, x FOR PEER REVIEW
M |F net −F| and M |F pro −F| increased with the increasing speed of the treadmill (p < 0.0 0.001, Table 2). The overall absolute mean difference increased by 23% from 3 to 5 the V2 technique, the overall CMC decreased by about 2% from 13 to 21 km/h. was independent of the incline of the treadmill (p = 0.042, Table 1). CMCpro increase 3 to 5° (p = 0.007, Table 1). The overall absolute difference increased by 33% from 1 km/h. M |F net −F| was dependent on the incline of the treadmill (p = 0.014, Table M |F pro −F| was independent of the incline of the treadmill (p = 0.577, Table 2).   On average, the M F net −F was lower than zero and the M F pro −F was greater than zero ( Figure 7) for both the DP and V2 techniques at any speeds and any inclines. The interaction effect (comparison * speed) was significant on the absolute mean force difference with the DP technique (p = 0.025) but not with the V2 technique (p = 0.165). In the DP technique, M |F net −F| was 24% lower than M |F pro −F| at 15 km/h (p = 0.028, Table 2). For the V2 technique, the overall M |F pro −F| was about 37% greater than M |F net −F| . The interaction effect (comparison * incline) was not significant on absolute mean force difference in the DP technique (p = 0.393) but was significant in the V2 technique (p = 0.016). In the DP technique, the overall M |F net −F| was about 39% lower than M |F pro −F| . With the V2 technique, M |F net −F| was significantly lower than M |F pro −F| at 3 • and 4 • (p ≤ 0.013, Table 2). Mean force difference over force producing phases in DP and V2 techniques (%BW). M F net −F represents the difference between F and F net and is calculated by F-F net . M F net −F lower than zero indicates that F net is greater than F. M F pro −F represents the difference between F and F pro and is calculated by F-F pro . M F pro −F greater than zero indicates that F pro is lower than F. With the DP technique, CMC pro was independent from the speed (p = 0.371, Table 1). However, CMC net decreased significantly at 21 km/h when compared to CMC net at 13, 15, and 17 km/h (p ≤ 0.046). The overall CMC increased by about 2% from 3 to 5 • . Both M |F net −F| and M |F pro −F| increased with the increasing speed of the treadmill (p < 0.001, p < 0.001, Table 2). The overall absolute mean difference increased by 23% from 3 to 5 • . With the V2 technique, the overall CMC decreased by about 2% from 13 to 21 km/h. CMC net was independent of the incline of the treadmill (p = 0.042, Table 1). CMC pro increased from 3 to 5 • (p = 0.007, Table 1). The overall absolute difference increased by 33% from 13 to 21 km/h. M |F net −F| was dependent on the incline of the treadmill (p = 0.014, Table 2), and M |F pro −F| was independent of the incline of the treadmill (p = 0.577, Table 2).

Discussion
The results of this study support our first hypothesis that the force-time curves of F net and F pro all give comparable shape with F in both techniques. In the DP technique, CMC pro ranged from 0.907 to 0.936, CMC net ranged from 0.883 to 0.955 and did not differ from CMC pro (Table 1). In the V2 technique, CMC pro ranged from 0.832 to 0.900, and CMC net ranged from 0.879 to 0.922 ( Table 1). The CMC depicting the similarity between waveforms and CMC close to 1 indicated that the curves involved were similar [21,22]. Therefore, the shapes of force-time curves of F pro and F net all showed similar to force-time curves of F, and both could be used to describe the shape of F during the poling phase of the DP technique and the kicking phase of the V2 technique. In addition, in the V2 technique, CMC net was 5% higher than CMC pro while changing the speed ( Table 1), indicating that the force-time curves of F net was more comparable to the force-time curves of F than F pro while using the V2 technique. Consequently, F net appears to be more appropriate for determining the trend of the forward acceleration in the V2 technique.
The results of this study partly support our second hypothesis that F net would give a considerable overestimation and F pro would be more accurate than F net when estimating F in both the DP and V2 techniques. In this present study, the F net had a considerable overestimation when estimating F in both the DP and V2 techniques, but F pro was not more accurate than F net in both techniques. The mean force differences over force curves between F pro and F (M F pro −F ), as well as F net and F (M F net −F ), were computed (Figure 7).
The mean force differences in this study indicated that, on average, F net would overestimate (M F net −F < 0, Figure 7) the F in both the DP and V2 techniques. F net was calculated from the GRF directly. The costs associated with the transformations of energy [24] between each segment and the elastic potential energy of the muscle were not taken out from F net . Thus, a considerable difference in F net and F may exist. The F pro underestimate (M F pro −F > 0, Figure 7) the F, but it was not more accurate than F net in both techniques. F pro was calculated by combining the GRF and the position of COM. The resultant GRF was subdivided into a translational component, which acted through the COM, and a rotational component, which was always perpendicular to the translational component [6,25]. Because the rotational component will not have a translational effect on the COM, when F pro was calculated, the rotational component was not involved. Therefore, the forward component of the translational component might underestimate the forward acceleration in both the DP and V2 techniques. The absolute mean force differences ( Table 2) were computed to evaluate which force-time curve was closer to the reference one. A smaller absolute mean force difference indicates a force-time curve closer to the reference curve and further shows a relatively higher accuracy. The results of this study showed that with both the DP and V2 techniques, the absolute mean force difference between F pro and F were greater than or have no difference with the absolute mean force difference between F net and F. This indicates that the force-time curves of F net were closer to or have no difference with the force-time curves of F. Thus, F pro was not more accurate than F net .
The results of this study do not support our third hypothesis that the approaches to calculate the F net and F pro would not be affected by the speed and incline of the treadmill in both techniques. The approaches to calculate the F net and F pro were all influenced by the speed or the incline of the treadmill. As there was a balance of forces under laboratory conditions with no air resistance, constant friction coefficient, and constant gravitational force, the total external force remained constant when the speed was changed. The gravity component parallel to the treadmill surface increased with the incline [19]; thus, more forces were needed at a steeper incline. It is impossible to have an exact reproduction of the reference value F; however, if the methods for calculating F net and F pro were independent from the speed and the incline of the treadmill, the CMC net , CMC pro , and the absolute mean force difference over the force-generating cycle should remain constant in both techniques. The CMCs in this study were somehow affected by the speed and incline of the treadmill in both the DP and V2 techniques. In addition, the results of this study showed that the absolute mean force differences between F net and F (M |F net −F| ) and between F pro and F (M |F pro −F| ) were all affected by the speed of the treadmill regardless of whether the DP or V2 technique was performed ( Table 2). The absolute mean force differences increased with increasing speed (Table 2), which means that although F pro and F net can be used to estimate the force-time curve of F, they do not remain stable when the speed changes. Thus, when investigating how F adapts to increasing speed by using F net or F pro , the increasing mean force differences should be considered. Both M |F net −F| and M |F pro −F| increased when the DP technique was used while increasing the incline of the treadmill (Table 2). However, when the V2 technique was used, M |F net −F| was affected by the increasing incline, and the significant increase was only found at the steepest incline (Table 2), but M |F pro −F| was not influenced by the incline of the treadmill. Thus, compared to F net , F pro was more stable when estimating F while changing the incline of the treadmill.
Therefore, when considering the whole poling phase in the DP technique, both F pro and F net are appropriate for estimating the trend of F. The similarity between the F pro and F is stable while changing the speed in the DP technique. However, F net has better accuracy than F pro when the speed and the incline is changed. When considering the whole kicking phase in the V2 technique, the trend of F net fits F better. However, the similarity between the F pro and F is stable in the V2 technique when the incline is changed. As the result in the DP technique, F net also has better accuracy than F pro in the V2 technique. There are some limitations of this study. The calculation of the COM is dependent on the assumed mass distributions. Although this has been proved to be suitable for estimating the position of COM in sports, it can still cause the golden standard of the reference to be inaccurate. In addition, the PFAs of the leg force and pole force were estimated points, and this may also have some effects on the accuracy of F pro . Furthermore, the added measurement equipment could have affected skiing performance.

Conclusions
The present study evaluated two approaches for estimating the total propulsive force on skier's COM. Both approaches can estimate the trend of the force-time curve of the propulsive force properly. Although both had a considerable overestimation; an approach by calculating the forward-directed horizontal component of 3D GRF is a more appropriate method due to a better accuracy. Future studies could investigate the contribution of skis and poles to forward COM acceleration by calculating the propulsive force from skis and poles separately. Moreover, as for the gliding phase that exists in XC skiing, the velocity at the end of the force generating phase is important. Future studies could also investigate the contributions of skis and poles to velocity change separately. Informed Consent Statement: Informed consent was obtained from all subjects involved in the study. Written informed consent was obtained from the patient(s) to publish this paper.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available because the data also form part of ongoing studies.

Acknowledgments:
The author Shuang Zhao wants to thank Christina Mishica and Ritva Mikkonen for language optimization. The authors would like to express their appreciation and thanks to the staff, athletes, and coaches for their participation, enthusiasm, and cooperation in this study.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Measurement and Calculation of Resistance Friction Coefficient of Roller Ski
The resistance friction coefficient of roller ski was measured on the treadmill surface using a custom-made friction measurement device (University of Jyväskylä, Jyväskylä, Finland, Figure A1) and calculated with the LabVIEW software package (National Instruments, Austin, TX, USA) before the measurement. A commercial force sensor (Raute precision TB5, Nastola, Finland) that measures the anterior-posterior force along the roller ski (F Y ) was contained in the friction measurement device. The friction coefficient between the treadmill surface and the roller ski was obtained by µ = F Y F Z , where F Z is the vertical force that equals the weight of the weight plate placed on the roller ski.
treadmill surface and the roller ski was obtained by μ = F Y F Z , where FZ is the verti that equals the weight of the weight plate placed on the roller ski. Figure A1. Custom-made friction measurement device.

Appendix B. Mechanical Principle of Translational Force
The motion of a rigid body under external forces can be reduced to (i) the acce of the COM and to (ii) the angular acceleration of the object around its COM. Thes the rate of momentum and angular momentum, respectively. Correspondingly space, the motion problem involves six degrees of freedom (DoF). These DoFs ca pressed with three components of a translational force and another three compo a moment. The forces and torques make a rigid body translate and rotate. It is wo ing that it is a modeler's decision to express the six DoFs with translational fo torques with respect to the COM. Accordingly, this is not the only, but rather a p option to model motion. The decomposition of an external force into translation and torque components acting on a rigid object is illustrated in Figure B1. An force F resultant ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ acts on point a of a rigid sphere. F resultant ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ is decomposed to the tran component F translational ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ that acts in the direction of the line joining the COM and The (displacement) vector from COM to point a is denoted by l. The rotational nent F rotational ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ of the force is perpendicular to vector l such that condition F r ⃗⃗⃗⃗ F rotational ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ + F translational ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ holds ( Figure B1a). Precisely, the same situation is expr terms of a translational force F translational ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ and torque τ, which is the product F rotational ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ( Figure B1b). Figure B1. Diagram of the mechanical principle of translational force.
As the principles of mechanics do not depend on the object-that is, Newto apply to all objects-the basic setting does not change from that of a single rigi However, in the case of joined bodies, the COM cannot be specified a priori, as it on the position of the parts in relation to each other (see Figure B2). When extern F resultant ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ act on the object (Figure B2a), the whole object translates, and each par

Appendix B. Mechanical Principle of Translational Force
The motion of a rigid body under external forces can be reduced to (i) the acceleration of the COM and to (ii) the angular acceleration of the object around its COM. These change the rate of momentum and angular momentum, respectively. Correspondingly, in a 3D space, the motion problem involves six degrees of freedom (DoF). These DoFs can be expressed with three components of a translational force and another three components of a moment. The forces and torques make a rigid body translate and rotate. It is worth noting that it is a modeler's decision to express the six DoFs with translational forces and torques with respect to the COM. Accordingly, this is not the only, but rather a practical, option to model motion. The decomposition of an external force into translational force and torque components acting on a rigid object is illustrated in Figure B1

Appendix B. Mechanical Principle of Translational Force
The motion of a rigid body under external forces can be reduced to (i) the acce of the COM and to (ii) the angular acceleration of the object around its COM. These the rate of momentum and angular momentum, respectively. Correspondingly, space, the motion problem involves six degrees of freedom (DoF). These DoFs ca pressed with three components of a translational force and another three compo a moment. The forces and torques make a rigid body translate and rotate. It is wo ing that it is a modeler's decision to express the six DoFs with translational for torques with respect to the COM. Accordingly, this is not the only, but rather a p option to model motion. The decomposition of an external force into translation and torque components acting on a rigid object is illustrated in Figure B1. An force F resultant ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ acts on point a of a rigid sphere. F resultant ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ is decomposed to the tran component F translational ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ that acts in the direction of the line joining the COM and The (displacement) vector from COM to point a is denoted by l. The rotational nent F rotational ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ of the force is perpendicular to vector l such that condition F r ⃗⃗⃗⃗ F rotational ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ + F translational ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ holds ( Figure B1a). Precisely, the same situation is expr terms of a translational force F translational ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ and torque τ, which is the product F rotational ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ( Figure B1b). Figure B1. Diagram of the mechanical principle of translational force.
As the principles of mechanics do not depend on the object-that is, Newto apply to all objects-the basic setting does not change from that of a single rigi However, in the case of joined bodies, the COM cannot be specified a priori, as it on the position of the parts in relation to each other (see Figure B2). When extern F resultant ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ act on the object (Figure B2a), the whole object translates, and each par with respect to each other. Because of this, the COM moves with respect to th However, the motion of the object still fulfils Newton's laws ( Figure B2b). Conse Figure B1. Diagram of the mechanical principle of translational force.
As the principles of mechanics do not depend on the object-that is, Newton's laws apply to all objects-the basic setting does not change from that of a single rigid object. However, in the case of joined bodies, the COM cannot be specified a priori, as it depends on the position of the parts in relation to each other (see Figure B2). When external forces − −−−− → F resultant act on the object (Figure B2a), the whole object translates, and each part moves with respect to each other. Because of this, the COM moves with respect to the parts. However, the motion of the object still fulfils Newton's laws ( Figure B2b). Consequently, nothing prevents one from decomposing the external forces into components of translational forces and torques with respect to the COM.