Propulsive Power in Cross-Country Skiing: Application and Limitations of a Novel Wearable Sensor-Based Method During Roller Skiing

Cross-country skiing is an endurance sport that requires extremely high maximal aerobic power. Due to downhill sections where the athletes can recover, skiers must also have the ability to perform repeated efforts where metabolic power substantially exceeds maximal aerobic power. Since the duration of these supra-aerobic efforts is often in the order of seconds, heart rate, and pulmonary VO₂ do not adequately reflect instantaneous metabolic power. Propulsive power (P_prop) is an alternative parameter that can be used to estimate metabolic power, but the validity of such calculations during cross-country skiing has rarely been addressed. The aim of this study was therefore twofold: to develop a procedure using small non-intrusive sensors attached to the athlete for estimating P_prop during roller-skiing and to evaluate its limits; and (2) to utilize this procedure to determine the P_prop generated by high-level skiers during a simulated distance race. Eight elite male cross-country skiers simulated a 15 km individual distance race on roller skis using ski skating techniques on a course (13.5 km) similar to World Cup skiing courses. P_prop was calculated using a combination of standalone and differential GNSS measurements and inertial measurement units. The method's measurement error was assessed using a Monte Carlo simulation, sampling from the most relevant sources of error. P_prop decreased approximately linearly with skiing speed and acceleration, and was approximated by the equation Pprop(v,v˙) = −0.54·v −0.71·v˙ + 7.26 W·kg⁻¹. P_prop was typically zero for skiing speeds >9 m·s⁻¹, because the athletes transitioned to the tuck position. Peak P_prop was 8.35 ± 0.63 W·kg⁻¹ and was typically attained during the final lap in the last major ascent, while average P_prop throughout the race was 3.35 ± 0.23 W·kg⁻¹. The measurement error of P_prop increased with skiing speed, from 0.09 W·kg⁻¹ at 2.0 m·s⁻¹ to 0.58 W·kg⁻¹ at 9.0 m·s⁻¹. In summary, this study is the first to provide continuous measurements of P_prop for distance skiing, as well as the first to quantify the measurement error during roller skiing using the power balance principle. Therefore, these results provide novel insight into the pacing strategies employed by high-level skiers.

Introduction

Current cross-country ski races can be categorized into two main events; sprint skiing (<1.8 km) and distance skiing (>10 and 15 km, female and males respectively). Furthermore, races can be held using free technique, where athletes typically choose to use ski skating techniques, or as classic technique races, where skiers are restricted to specific sub-techniques (herringbone, diagonal stride, double poling, kick double poling). Due to these restrictions, the average race speed in free technique events is typically about 10% higher than in classical technique races (Bolger et al., 2015). Regardless of events, the race course regulations specify that courses should contain approximately equal parts of uphill, downhill, and flat terrain, to “test the skier in a technical, tactical and physical manner” (FIS Cross-Country Homologation Manual, June 2017).

Regardless of race distance, cross-country skiing is an endurance sport that demands an exceptionally high aerobic energy turnover, in addition to high movement efficacy. This is underlined by the fact that elite cross-country skiers have among the highest maximal oxygen consumptions of any sports (Sandbakk and Holmberg, 2014; Haugen et al., 2017), typically ranging from 80 to 90 and 70 to 80 mL·kg⁻¹·min⁻¹ for world class males and females, respectively (Ingjer, 1991; Losnegard and Hallén, 2014; Sandbakk et al., 2016a). In addition, several studies also indicate substantial anaerobic turnover rates during a race, a phenomenon attributed to the large variations in course inclination. Moreover, skiers typically choose to increase their metabolic power in uphill terrain (Karlsson et al., 2018), often attaining a metabolic power that substantially exceeds their peak aerobic power (estimated at 110–160% of VO_2peak Norman and Komi, 1987; Sandbakk et al., 2011; Karlsson et al., 2018). These repeated supra-aerobic efforts vary in duration from seconds to minutes, and incur an oxygen debt that must be recovered in the downhill or flat sections (Sandbakk and Holmberg, 2014; Karlsson et al., 2018). Such transient changes in energetic demand are not well reflected by measurements of pulmonary VO₂ or heart rate, because both have a blunted response due to the use of local oxygen stores and anaerobic energy pathways. Hence, both parameters behave as if they passed through a lagged low pass filter and remain high (85–95% of their peak values) throughout the race (Welde et al., 2003; Bolger et al., 2015; Karlsson et al., 2018).

The combination of high and sustained aerobic energy turnover with repeated supra-aerobic efforts distinguishes cross-country skiing from many other endurance racing sports, where the work rate is relatively constant and requires measurement of parameters that reflect the instantaneous energy demands in a competition setting. A frequently used parameter that often corresponds well with instantaneous energy requirement is the propulsive power (P_prop) generated by the athlete. For some endurance sports, like cycling, P_prop can be measured directly, and metabolic energy requirements are approximately linearly related to P_prop (Ettema and Lorås, 2009). Hence, if P_prop can be linked to metabolic power in skiing (Millet et al., 2003; Sandbakk et al., 2011; Karlsson et al., 2018), in-field measurements of P_prop would be useful to further our understanding of the physiological demands experienced by a competitive skier. Therefore, several studies have attempted to calculate P_prop during cross-country skiing, based either on position measurements (Sandbakk et al., 2011; Swarén and Eriksson, 2017), or simulation of skiing performance (Moxnes and Hausken, 2008; Carlsson et al., 2011; Moxnes et al., 2013, 2014; Sundström et al., 2013). These studies all use the principle of power balance, as outlined by van Ingen Schenau and Cavanagh (1990). However, no studies are available where P_prop has been measured continuously throughout a cross-country ski race with a duration longer than sprint skiing. Furthermore, no previous studies have critically evaluated the accuracy achieved when applying the power balance principle to cross-country skiing. The aims of this study were (1) to develop a procedure for estimating the propulsive power generated during roller-skiing using small non-intrusive sensors (GNSS and IMUs) attached to the athlete and evaluate its limitations; and (2) to utilize this procedure to determine the propulsive power generated by high-level skiers during a simulated distance race.

Theoretical Background

As stated by van Ingen Schenau and Cavanagh (1990), the P_prop is equal to the rate of change in mechanical energy (E_mech) of the system and the work done on the environment (W_env). In cross-country skiing P_prop is customarily estimated by modeling the skier and his/her equipment as a point mass (Moxnes and Hausken, 2008; Carlsson et al., 2011; Moxnes et al., 2013, 2014; Sundström et al., 2013; Swarén and Eriksson, 2017). Under this assumption, the mechanical energy is the sum of translational kinetic energy and the gravitational potential energy. Work done on the environment is primarily due to ski/snow-friction forces (or rolling resistance, denoted F_f) and the aerodynamic drag force (F_d). This is summarized in Equation 1:

$\begin{array}{cc}\begin{array}{l}{P}_{\text{prop}}={\dot{E}}_{\text{mech}}+{\dot{\text{W}}}_{\text{env}}\\ =\frac{\left(m\dot{v}-{\mathbf{\text{F}}}_g-{\mathbf{\text{F}}}_{\text{d}}-{\mathbf{\text{F}}}_{\text{f}}\right)·\mathbf{\text{v}}.}{\text{Point mass assumption}}\end{array}& (1)\end{array}$

In Equation 1 m refers to the total mass of the system (the sum of body mass and equipment), v is the velocity of the center of mass (COM), $\dot{v}$ is the COM acceleration, and F_g is the gravitational force. Furthermore, the magnitude of the propulsive force (F_prop), i.e., the force in the skiing direction that is not due to gravity or frictional forces (air drag, ski/snow-friction or rolling friction) is calculated using F_prop = P_prop·|v|⁻¹ (Carlsson et al., 2011).

For skiing and roller skiing applications F_f has commonly been modeled using the Amonton-Coulomb equations (Carlsson et al., 2011; Moxnes et al., 2013, 2014; Sundström et al., 2013; Swarén and Eriksson, 2017). This friction model is attractive because of its simplicity, but it is unable to capture complex ski-snow interactions (Bowden and Hughes, 1939; Buhl et al., 2001; Theile et al., 2009), and F_f may change considerably over the course due to changing snow conditions. This challenge can be partially overcome using roller skis, which have a more constant coefficient of rolling resistance, except during the warm-up period (Ainegren et al., 2008). Another limitation is that during turns or when employing ski-skating techniques the skis (or roller skis) are actively pushed in the mediolateral direction, which causes shear forces. One study has investigated how the rolling resistance of a roller ski was affected by shear forces occurring during the ski push-off (Sandbakk et al., 2012). They concluded that the ratio of the rolling resistance force to the vector sum of shear and compression forces varied by <2% for ski angles up to 45°. This finding suggests that when evaluating rolling resistance during roller ski skating, shear forces (caused by ski push-off or centripetal forces during a turn) can be added to the normal force, at least for the wheel type assessed by Sandbakk et al. (2012).

Work done against the aerodynamic drag force has conventionally been modeled using the drag equation from fluid dynamics (Carlsson et al., 2011; Moxnes et al., 2013, 2014; Sundström et al., 2013; Swarén and Eriksson, 2017):

$\begin{array}{cc}{F}_d =-\frac{1}{2}\rho v_f^2·A{C}_D(\text{Re})·{\widehat{v}}_f.& (2)\end{array}$

In Equation 2, ρ denotes the air density, v_f denotes the velocity relative to the air, and $\widehat{\mathbf{\text{v}}}$_f the unit vector along v_f. The drag area (AC_D) is the product of the frontal area A of the athlete and equipment and the drag coefficient (C_D), which depends on the shape and surface material properties of the object in the air flow. C_D also depends upon the Reynolds number:

$\begin{array}{cc}\begin{array}{r}\begin{array}{c}\text{Re}=\frac{\rho L{v}_f}{\mu }.\end{array}\end{array}& (3)\end{array}$

In Equation 3, ρ denotes the fluid's density, L is the characteristic length of the object, and μ is the dynamic viscosity. C_D is relatively constant when the flow is turbulent, except for a sharp drop when the boundary layer transitions from semi-turbulent to fully turbulent flow, resulting in a narrower wake behind the object (Spurk and Aksel, 2008). This phenomenon usually occurs at Re around 2 × 10⁵ and is well studied for simple blunt bodies (Achenbach, 1968; Spurk and Aksel, 2008). The drop in C_D has also been shown to exist for athletes while roller skiing (Spring et al., 1988) and during wind tunnel simulations of ice skating (van Ingen Schenau et al., 1982), where the transitions occurred at about 5 and 10 m·s⁻¹, respectively. The magnitude of the change in C_D varies considerably between different studies. Data from Achenbach (1968) and van Ingen Schenau et al. (1982) showed that C_D was reduced to about 30–40% of its quasi-stable value for Re < 10⁵, while data from Spring et al. (1988) implied a decrease to about 10%. There are two more challenges when estimating F_d from Equation 2: (i) wind velocity must be known in order to find v_f; and (ii) the drag area depends on the skier's posture and (indirectly) on skiing speed, through the Reynolds number. Previous studies have estimated AC_D as constant (Swarén and Eriksson, 2017), or as a step or smooth function of skiing speed using allometric scaling based on body mass (Sundström et al., 2013; Moxnes et al., 2014), and only one study has investigated the effect of non-zero wind velocity (Moxnes et al., 2014). These simplifications are often necessary because direct measurement of instantaneous wind field and drag area are challenging or impossible in the field. Nonetheless, the error arising from these approximations has rarely been addressed.

Materials and Methods

Participants and Study Design

The data presented in this study were collected over three test days on the roller skiing course at Holmenkollen, Oslo, Norway. The topography of the course is similar to race courses used in competitive cross-country skiing on snow (height difference 51 m, maximum climb 32 m, total climb 166 m). Eight skiers (seven cross-country skiers; FIS point range of 13–117, and one biathlete) volunteered for the study and gave their written consent to participate. The study was approved by the ethics committee at the Norwegian School of Sport Sciences and the Norwegian Centre for Research Data, and was conducted according to the Declaration of Helsinki.

The participants were asked to complete a test race consisting of three laps of a 4.5 km course in the shortest time possible. The test race was arranged as a time trial and the participants were instructed to the use skating techniques. All participants used the same model of roller skis (Swenor Skate Long, wheel type 2), and wore tight-fitting clothing. Each participant was equipped with two identical position tracking devices (Catapult Optimeye S5, mass 67 g) consisting of a 10 Hz standalone GNSS-module and a 9-axis inertial measurement unit (accelerometer, gyroscope and magnetic field measurements). One of the receivers, carried in a tight-fitting vest, was positioned at approximately the level of the third thoracic vertebra, while the other was taped laterally on the thigh approximately 10 cm inferior to the trochanter (Figure 1). Both units were attached so that the inertial sensors' local coordinate systems (xyz) were approximately aligned with the x-axis directed the mediolaterally to the right (in the skiing direction) and the y-axis in the anterior direction. Prior to the test race the weight of the athletes and equipment (including ski boots, roller skis, ski poles and helmet) was recorded, and the participants performed calibration measurements for the drag area model, as described in section Frontal Area. After completing these measurements, the athletes warmed up for approximately 30 min before the race started.

FIGURE 1

Figure 1. Positioning of the two IMUs on the athletes' body and motivation for the frontal area model described in Equation 5. One IMU was positioned approximately at the level of the third thoracic vertebra, the other was taped laterally on the thigh approximately 10 cm inferior to the trochanter. The axes were aligned with the x-axis (blue) in the mediolateral direction and the y-axis (red) in the anterior direction. The z-axis (green) was aligned with gravity when the athletes were in a standing posture, as described in the text. The frontal areas of the torso and thighs scale approximately with the cosine of the pitch angles θ _k = atan(a_{y, k}/a_{z, k}), where a_{y, k} and a_{z, k} are the smoothed accelerometer outputs along the y and z-axis respectively. k represents the thigh or torso location.

To allow more accurate determination of vertical position, the GNSS positions from the receivers carried by the athletes were mapped onto a common trajectory that was determined from kinematic position tracking using a more accurate differential GNSS (receiver: Alpha-G3T, antenna: GrAnt-G3T, Javad, USA). These measurements had an expected accuracy <5 cm (Gilgien et al., 2014b) when double difference ambiguities were fixed. The ambiguities could not be solved for five short sections of the course, resulting in a substantially reduced accuracy in these sections (expected errors >1 m Gilgien et al., 2014b). These sections are clearly marked in the results.

Environmental Conditions

Air temperature and pressure were recorded during the race for each day of testing. Wind velocities (measured 10 m above the ground) were retrieved hourly from weather stations located approximately 2.3 km (Tryvann) and 3.7 km (Blindern) from the roller skiing course (www.met.eklima.no, Meteorological Institute of Norway, Oslo, Norway). The wind velocities were averaged and corrected to 1 m height above the ground using the log wind profile relationship (Oke, 1978) with a roughness length of 0.1 m, resulting in wind speeds exactly half of the data measured 10 m above the ground.

Data Processing and Filtering

The GNSS positions and inertial sensor data from the receivers mounted on the athletes were exported using Catapult Sprint software version 5.1.7 at sampling frequencies of 10 and 100 Hz, respectively. Differential GNSS positions for the reference trajectory were calculated using the kinematic algorithm of the geodetic post-processing software Justin (Javad, San Jose, CA, USA). All other analyses were performed using Matlab R2017a (MathWorks, Natick, MA, USA).

The reference trajectory was resampled to equidistantly (1 meter) spaced points, and then filtered using smoothing splines weighted by their fixed/float status (Skaloud and Limpach, 2003) using a smoothing parameter of p = 0.02. The GNSS positions obtained from the receivers mounted on the athletes were filtered with a second order bidirectional Butterworth low pass filter with a cutoff frequency of 0.3 Hz, which removed the high frequency components that could not be attributed to the center of mass trajectory due to the receiver's antenna location. The cutoff-frequency was chosen based on a frequency-analysis of motion capture data from skiing using the V1 and V2 techniques. Only the GNSS positions from the receiver on the torso were used for skier position calculation. Vertical position (z) was obtained by mapping the horizontal plane position from the standalone GNSS receiver carried by the athlete onto the 3D reference trajectory. This was achieved by minimizing horizontal plane Euclidean distance from the position measurement of the standalone GNSS receiver to any point on the reference trajectory. After filtering and projection of the standalone GNSS positions onto the reference trajectory the trajectories were down-sampled to a frequency of 1 Hz prior to work rate calculations, to limit the computational load. Finally, external work rate was filtered using a 5 second bidirectional moving average filter.

Mechanical Energy

For calculation of mechanical energy, a point mass m equal to the combined mass of the athlete and his equipment was utilized. With this approach the gravitational potential energy is mgz, where g = 9.81m·s⁻² is the acceleration due to gravity and z the vertical position. The kinematic energy is 0.5m|v|², where v is the skiing velocity. The skiing velocity was determined by differentiating the horizontal plane positions from the standalone GNSS receiver carried by the athlete and the vertical component of velocity was calculated using the using the mapped vertical position. Velocity was calculated using a five-point finite difference algorithm (Gilat and Subramaniam, 2008).

Rolling Resistance

Rolling resistance was measured individually for each athlete using a towing test on the roller skiing treadmill at the Norwegian School of Sport Sciences, as described by Hoffman et al. (1990). The coefficient of rolling resistance (C_rr) was 0.0225 ± 0.0009 (mean ± standard deviation). The rolling resistance of one of the pairs of roller skis on asphalt (determined by a towing test) was the same as on the treadmill surface, in agreement with the findings of Myklebust (2016), and the former were therefore used in the calculations of work against rolling resistance. The roller skis were assumed to move the same distance as the GNSS antenna, and centripetal forces (F_c, caused by the course's curvature) were added to the normal force opposing gravity (N_g), following the findings of Sandbakk et al. (2012). Hence, the work rate against rolling resistance was estimated using the following equation:

$\begin{array}{cc}\begin{array}{r}\begin{array}{c}{\dot{W}}_f=C_{rr}·\left|{\text{N}}_g+{\text{F}}_c\right|·v.\end{array}\end{array}& (4)\end{array}$

In Equation 4, N_g was defined as minus the component gravity perpendicular to the course's normal vector, and the course was assumed to be level in the mediolateral direction. Centripetal force was calculated using F_c = m|v|²·K, where K = −v×(v×$\dot{\text{v}}$)/v⁴ is the track curvature vector (Dooner, 2012).

Air Drag Model

Air drag was determined from Equation 2, with the drag area as defined in the next two paragraphs. Air density ρ_air = p·$R_{\text{specific}}^{-1}$T⁻¹ was calculated from measured air pressure using the ideal gas law with R_specific = 287.058 J·kg⁻¹·K⁻¹, and p and T equal to the measured air pressure and temperature, respectively. Dynamic viscosity (μ_air) was calculated as a function of air temperature using Sutherland's formula, as described by Canuto et al. (2007). Wind velocity was defined as described in section Statistics and was subtracted from the skiing velocity vector. Finally, the power dissipated through air drag was the dot product F_d·v.

Frontal Area

The skiers' frontal area A was approximated continuously during the entire trial using the accelerometer data provided by the sensors mounted on the torso and thigh. First, the accelerometers' coordinate frames were rotated so that the z-axis was parallel to the gravity vector during a static standing pose (Figure 2A, pose 1), and the x-axis was directed mediolaterally. This was achieved by performing two successive rotations, the first canceling lateral tilt and the second canceling forward tilt (Myklebust et al., 2015). Second, the signals were filtered with a 2 second bidirectional moving average filter. Third, the smoothed pitch angles θ _thigh and θ _torso were calculated using tan(θ) = a_y/a_z, where a_y and a_z denote the smoothed accelerometer outputs (in a standing posture) along the forward and vertical directions. With this definition, the pitch angle represents a rotation about the IMU's x-axis that will align the y-axis with the horizontal plane. As shown in Figure 1, the frontal areas of the torso or thigh segments will scale approximately with the cosine of this angle. Therefore, θ _thigh and θ _torso were used to predict the frontal area A of the athlete using the following equation:

$\begin{array}{cc}\begin{array}{r}\begin{array}{c}\frac{A-A_{\text{equip}}}{A_{0,j}}={\beta }_0+{\beta }_{1}cos({\theta }_{torso})+{\beta }_{2}cos({\theta }_{thigh}).\\ \end{array}\end{array}& (5)\end{array}$

Here A_{0, j} denotes the frontal area of athlete j while standing upright, as shown in Figure 2A, and A_equip = 0.045 m² was a constant that was added to represent the average frontal areas of the roller skis and ski poles. The three parameters β_{0, 1, 2} were determined using multiple linear regression with cos(θ _{thigh, (i, j)}) and cos(θ_{torso, (i, j)}) as predictors, and 56 frontal areas A_{(i, j)} (i.e., seven frontal areas per participant) with different postures as the dependent variable (Figures 2A,C). The 56 frontal areas were calculated from digital images (resolution 3,264 × 4,928 pixels) taken of the skiers prior to the trial, and the pixel size was determined using an object of known length placed directly lateral to the athlete. The characteristic length L_j in Equation 3 was defined as the width of athlete j in the pelvis/abdomen region and was determined from the first of the 7 images (Figure 2A).

FIGURE 2

Figure 2. (A) Postures used to calculate frontal area. The first pose for each athlete was used to calculate the characteristic length L (Equation 3), which was defined as the width at the height of the pose's center of mass (green line). (B) Model for the frontal area A based on the postures in (A) and Equation 5. The data points are enumerated to match the postures in (A), and show the measurements from one subject. Red dots indicate measurements from all other subjects. (C) Model (Equations 7, 8) for the drag coefficient of a skier as a function of the Reynolds number (solid black line). Blue squares: data from Achenbach (1968); Green diamonds: data from van Ingen Schenau et al. (1982); Red circles: data from Spring et al. (1988); Dashed blue line: model fitted to Achenbach's measurements. The upper and lower x-axis were related by Equation 3 using atmospheric conditions as specified in van Ingen Schenau et al. (1982) and L = 0.3 m. Shaded areas indicate the models' standard deviations, and were obtained using Monte Carlo simulations (N = 5,000) sampled from the distributions specified in the statistics section.

Determination of the Frontal Area by Allometric Scaling

We also tested a drag area model that simplifies data collection and analysis, because it does not require measurements of frontal areas or measurements from inertial sensors. First, the estimated frontal areas determined from Equation 5 were normalized by m^2/3, where m was body mass and 2/3 is the allometric scaling exponent (Günther, 1975; Bergh, 1987). Second, the normalized frontal areas were assumed to be a sigmoid function of skiing speed v, as previously assumed by Moxnes et al. (2014) and Sundström et al. (2013). To model this sigmoid shape, we defined a logistic function:

$\begin{array}{cc}\begin{array}{r}\begin{array}{c}\frac{A(v)}{m^{2/3}}={\gamma }_1-\frac{{\gamma }_2}{1+e^{-\left(\frac{v}{{\gamma }_4}-{\gamma }_3\right)}}.\\ \end{array}\end{array}& (6)\end{array}$

The four-parameter vector γ was determined by minimizing the sum of the squared residuals from the normalized estimates of drag area (A·m^−2/3) determined from Equation 5 using the Levenberg-Marquardt algorithm. Hence, once γ was established, the only necessary input was body mass and skiing speed.

Drag Coefficient

The drag coefficient C_D was modeled as a function of the Reynolds number (Re) using a logistic function fitted to the data from Achenbach (1968). Specifically, the four-parameter vector α was determined by minimizing the sum of the squared residuals between measurements and the model in Equation 7 (below) using the Levenberg-Marquardt algorithm (implemented in Matlab's Statistics and Machine Learning toolbox):

$\begin{array}{cc}\begin{array}{r}\begin{array}{c}{C_D}^{\prime }(\text{Re})={\alpha }_1-\frac{{\alpha }_2}{1+e^{-\left(\frac{\text{Re}}{{\alpha }_4}-{\alpha }_3\right)}}.\\ \end{array}\end{array}& (7)\end{array}$

Only measurements with Reynolds numbers in the range [4 × 10⁴, 5 × 10⁵], which corresponds to conditions relevant to cross-country skiing (speeds in the range 2–25 m·s⁻¹, assuming μ·ρ⁻¹ = 1.5 × 10⁻⁵m²·s⁻¹ and L = 0.30 m) were used for the fitting procedure (Figure 2C). In the second step, the model in Equation 7 was scaled to match the mean drag coefficients of speed skaters measured at 12 m·s⁻¹ by van Ingen Schenau et al. (1982). This was done by calculating the Reynolds number (Re_Schenau) for the conditions described in van Ingen Schenau et al. (1982), and then applying the following equation:

$\begin{array}{cc}\begin{array}{r}\begin{array}{c}{C}_D(\text{Re})=\frac{C_{D,\text{Schenau}}}{{C_D}^{\prime }({\text{Re}}_{\text{Schenau}})}{C}_D^{{}^{\prime }}(\text{Re}).\end{array}\end{array}& (8)\end{array}$

Tuck Position

When in the tucked position, skiers do not generate any significant propulsive force. As all cross-country skiing techniques (except the tucked position) require substantial rotation of both the thorax and thigh segments, measurements of the segments' angular rates were used to determine when the athletes were in the tucked position. Specifically, the squared magnitude of the angular velocity vector from the gyroscopes was used as a decision criterion. These signals were filtered with a 2 second bidirectional moving average filter, and athletes were defined to be in the tuck position when both devices (on thigh and thorax) showed values smaller than 5,000^°2·s^−2. This threshold was determined by inspecting the signal distributions (Figure 3). P_prop was set to zero when the athletes were in the tucked position.

FIGURE 3

Figure 3. Histogram of the squared magnitude of gyroscope measurements (angular velocity) during skiing for all athletes. The distinct distribution minimum at about 5,000^°2s⁻² indicates the transition from the tucked position, where movements of both torso and thigh are small, to skiing techniques that generate propulsion. In the current study, the athletes were assumed to be in the tucked position when |ω|² of both torso and thigh were <5,000^°2·s⁻².

Statistics

Confidence intervals for the calculated work rates were determined using a Monte Carlo approach. The effect of changing wind velocity and errors in drag area and rolling resistance were modeled by sampling from distributions as described in this section. 2,500 Monte Carlo samples were used for each of the participants.

Wind speed was assumed to be Rayleigh distributed (Dorvlo, 2002) with an expectation value equal to the hourly average wind speed after correction for height above ground, and averaged over the two measurement stations. Wind direction was assumed to be normally distributed with an expectation value equal to the hourly average and a standard deviation of 25 degrees.

Errors in frontal area estimates (Equation 5) were assumed to be ~N(0, ${\sigma }_{\text{resid}}^2$), σ_resid being the standard deviation of the model residuals. The distribution of the coefficients α in Equation 7 was assumed to be ~N(α_opt, ${\sigma }_{\text{opt}}^2$), where α_opt and σ_opt were the coefficient estimates and covariance matrix obtained from the optimization procedure. C_{D, Schenau} was defined to be ~N(0.872, 0.079²), which corresponds to the mean and standard deviation of the findings in Table 1 of van Ingen Schenau et al. (1982).

The coefficient of rolling resistance was assumed to be normally distributed with an expectation value equal to the measured values from the treadmill towing test and standard deviation 2.3 × 10⁻³, which was the standard error of the measurements from the towing test on asphalt. Asphalt values were used for the standard deviation, since these were assumed to better represent the variability of field conditions.

The method's accuracy was defined as the pooled standard deviation of P_prop and F_prop obtained from the Monte Carlo simulation. The method's sensitivity was assessed by comparing inter-athlete and intra-athlete differences in F_prop to the method's accuracy. This was achieved by calculating the empirical cumulative distribution function of the inter-athlete or intra-athlete standard deviation (i.e., the lap-to-lap variability) of F_prop. The method was considered suitable to discriminate differences in F_prop if the standard deviation was greater than the measurement accuracy for >90% of the measurements. To test the method's sensitivity in optimal conditions, the Monte Carlo simulation was run again with zero wind speed. All other parameters were kept constant.

Results

Physiological Aspects

A graphical presentation of P_prop as a function of distance traveled is provided in Figure 4, which clearly shows that there were substantial variations in P_prop throughout the test race. P_prop was 3.35 ± 0.23 W·kg⁻¹ when averaged over the race duration, and 4.18 ± 0.41 W·kg⁻¹ when omitting measurements where the athletes were in the tucked position (termed active propulsive power, Swarén and Eriksson (2017)). Peak P_prop was typically attained on the last major ascent during the final lap (at 3,600–3,730 m from the start of the lap in Figure 4). The average of the athletes' peak P_prop was 8.35 ± 0.63 W·kg⁻¹. When comparing P_prop at the same positions along the course, lap-to-lap differences were small, except for a distinct starting spurt for the initial 200 meters of the first lap (Figures 4B,C). P_prop was also higher in the last two uphills of the third lap (3,600–3,730 and 4,120–4,170 m) compared to the first lap, indicating an end spurt (Figure 4C).

FIGURE 4

Figure 4. (A) Propulsive power normalized to body mass plotted over distance along the course for the three laps (lap 1 red, lap 2 green, lap 3 blue). Colored regions around the solid lines show the standard deviations from the Monte Carlo simulation. Vertical gray shading indicates regions where double differenced ambiguities were float. The Monte Carlo simulation does not account for the reduced accuracy in these regions. Negative values of P_prop occurs when F_prop < 0 and could be caused by either active breaking of the athlete, or by measurement error. (B) Mean difference in propulsive power between lap 2 and lap 1. Colored shaded region shows the 95% CI. (C) Same as B, but for the difference between lap 3 and lap 1. (D) Altitude profile of the competition course. Black regions correspond to the regions where double differenced ambiguities were float.

Overall, P_prop showed an approximately linear relationship with skiing speed (Figure 5A), except for speeds where the skiers were in the tucked position. The transition into the tucked position occurred at skiing speeds ~9 m·s⁻¹ (Figure 6). Furthermore, P_prop was dependent on the acceleration along the skiing direction; a positive acceleration correlated with smaller P_prop, and negative accelerations with higher P_prop. This is consistent with observations from Figure 4A, where P_prop appeared to decline slightly from its initial values in the longer ascents (450–600, 1,200–1,330, 2,030–2,450 m from the start, Figure 4D). This implies that P_prop was higher when the skier decelerated at the start of a climb, and lower when accelerating over the top of a climb. A multiple least squares regression fit on P_prop using skiing speed and acceleration as predictors had the coefficients −0.54 N·kg⁻¹, −0.71 N·s·kg⁻¹ and intercept 7.26 W·kg⁻¹ (Figure 5A), and a coefficient of determination R² = 0.403. This indicates that a substantial fraction of the variability in P_prop could not be explained by skiing speed and acceleration alone. In Figure 5B skiing speed is plotted against the course inclination. The line predicting steady state skiing speed based on the least squares fit P_prop(v,$\dot{v}$) was also added to the figure (see figure legend for details). The data points where the skier was close to a steady state skiing speed (i.e., |$\dot{v}$| is small) fell along the line suggested by the least squares fit. Data points far from the line suggested by the least squares fit was typically from periods with large positive acceleration (points below the line, colored red in Figure 5B) or negative acceleration (points above the line, colored blue in Figure 5B).

FIGURE 5

Figure 5. (A) Propulsive power (normalized to body mass) was approximately linearly related to skiing speed (v) and acceleration ($\dot{v}$) in the skiing direction. Data from all athletes are included in the figure, and data points are color coded by $\dot{v}$, as indicated by the colorbar above the figure. The magenta-colored line is the least squares regression fit to samples where the athlete was not in the tucked position. The shaded region indicates the measurement error (SD), see Figure 7 for details. (B) Plot of skiing speed vs. the course inclination. Data points are color coded as in panel (A). The black lines indicate constant propulsive power (in W·kg⁻¹), assuming constant skiing speed ($\dot{v}$ = 0). The magenta line shows the steady state skiing speed obtained from the regression line in (A). It was defined as the real root of the 3rd degree polynomial obtained by replacing P_prop with P_prop(v,$\dot{v}$) in Equation 1, and assuming a constant drag area of 0.55 m² and average body mass 77.1 kg. For v > 9 m·s⁻¹ the line was defined to follow the zero-P_prop iso-line. The figure clearly shows that data points where $\dot{v}$ ≈ 0 (green color) are distributed close to the steady-state speed line, as expected. Data points below the steady-state line have $\dot{v}$ > 0 (red), and data points above the line have $\dot{v}$ < 0. This is mainly attributed to the athletes' inertias.

FIGURE 6

Figure 6. (A) Frontal areas (Â) estimated from the accelerometer outputs plotted over skiing speed, normalized by body mass (using an allometric scaling coefficient of 2/3). Red dots represent measurements where the skiers were in the tucked position. Blue line: logistic function fitted to the data (Equation 6). Athletes typically assumed the tucked position at skiing speeds >9.1 m·s⁻¹. (B) Drag areas normalized to body mass (·C_D) found in the current study plotted over skiing speed for all athletes. The blue line shows the product of the frontal area model from Equation 6 (plotted in panel A) and the drag coefficient (Equation 8, Figure 2). Models used by Sundström et al. (2013), Moxnes et al. (2014), and Swarén and Eriksson (2017) are included for comparison.

Methodological Aspects

In Figure 6, Frontal areas (normalized to body mass) are plotted vs. skiing speed. The frontal area changed from the typical standing posture for v < 8 m·s⁻¹ to a tucked-position for v > 10 m·s⁻¹ (Figure 6A). Due to the dependency of C_D on the Reynolds number, the behavior was more complex for the drag area, which continuously decreased almost throughout the range of skiing speeds observed during the race (Figure 6B).

The Monte Carlo simulation showed that errors in F_prop were on average 5.1 · 10⁻² N·kg⁻¹ and increased slowly with skiing speed (Figure 7A). The least squares regression line of F_prop with respect to skiing speed had a slope of 2.8 ·10⁻³ N·s·m⁻¹·kg⁻¹ and intercept 3.9 · 10⁻² N·kg⁻¹. Therefore, the error in P_prop increased curvilinearly with skiing speed, approximately following the expression 2.8·10⁻³·v² + 3.9· 10⁻²·v (Figure 7B). Hence, estimates of P_prop were most accurate at low skiing speeds (~0.09 W·kg⁻¹ at 2.0 m·s⁻¹, close to the lowest measured speeds) and less accurate at high speeds (~0.58 W·kg⁻¹ at 9.0 m·${\text{s}}_{,}^{-1}$ close to the speed when most athletes transitioned to the tucked position).

FIGURE 7

Figure 7. (A) Standard deviation of propulsive force (normalized to body mass) calculated from the Monte Carlo simulation plotted over skiing speed. The red dots are measurements when the athletes were in the tucked position. In a later step propulsive power (and force) was defined to be zero for these samples, but they are included in this figure to illustrative purposes. The blue line is the least squares linear regression line fitted to the data points where the athletes were not in the tucked position. (B) same as (A), but for P_prop. Since P_prop = F_prop·v, the regression line from (A) multiplied with v fits well to the data.

The method's applicability to discriminate between instantaneous inter-athlete and intra-athlete differences in F_prop is shown in Figure 8. The results show that with the test conditions in the current study, neither inter-athlete or intra-athlete differences could be reliably detected. This conclusion holds true even in zero-wind conditions, where 29.9% (intra-athlete) and 79.9% (inter-athlete) of the measurements contained differences that exceeded the measurement accuracy.

FIGURE 8

Figure 8. In order to discriminate the F_prop or P_prop generated by high-level athletes, the methods accuracy must be better than typical athlete-to-athlete or within-athlete differences. To that end, this figure displays the empirical complementary cumulative distribution function (ECCDF) of typical athlete-to-athlete differences in F_prop measured at the same position along the course (blue line). More specifically, it shows the between-athlete standard deviation of F_prop evaluated at every integer meter along the course, excluding measurements where the athletes were in the tucked position. The red line is the ECCDF of the within-athlete lap-to-lap standard deviation. The solid vertical line at 0.0568 N·kg⁻¹ shows the mean standard deviation of F_prop from the Monte Carlo simulation, which reflects the methods typical accuracy under the specified conditions. For comparison, the dotted vertical line at 0.0385 N·kg⁻¹ shows the mean standard deviation of F_prop from another Monte Carlo simulation assuming zero-wind conditions. With the measured wind conditions, inter-athlete SD was greater than the typical measurement accuracy only for 39.1% of the course, and intra-athlete SD only for 11.9% of the course. Under the zero-wind assumption, the inter- and intra-athlete differences were greater than the measurement accuracy for 79.9 and 29.9% of the course, respectively.

For all results reported above, P_prop was calculated using vertical position measurements obtained by mapping the positions onto a reference trajectory of dGNSS measurements and using frontal area estimates from two accelerometers positioned at the thorax and thigh. Omitting the mapping procedure (by using the vertical position measurements of the standalone receivers carried by the athletes) yielded a root mean square (RMS) deviation of 1.25 W·kg⁻¹ from the more accurate method of using vertical position mapped onto the dGNSS measurements. This is about 30% of mean active propulsive power. Example data for these calculations from one of the athletes are shown in Figure 9. Omitting the accelerometer measurements by using frontal areas obtained from allometric scaling (Equation 6, Figure 6) yielded an RMS deviation of 0.18 W·kg⁻¹, or about 4% of mean active propulsive power (Figure 9).

FIGURE 9

Figure 9. (A) Comparison of propulsive power calculations using mapping on dGNSS reference (blue line), calculations using only measurements from the standalone GNSS receiver (green line), and using the simplified drag area model based on body mass and skiing speed (Equation 6, red line). Data are from one lap for a single athlete. Vertical gray shading indicates regions where double differenced ambiguities were float. (B) Altitude profile of the competition course. Black regions correspond to the regions where double difference ambiguities were float.

Regression and curve fitting parameters for the frontal area model (Equation 5), drag coefficient model (Equation 8) and allometric scaling of frontal areas (Equation 6) are presented in Table 1.

TABLE 1

Table 1. Results for the regression parameters (β) for frontal area calculation (Equation 5), parameters (α) for the drag coefficient model (Equation 7), and parameters (γ) for prediction of frontal areas using body mass and skiing speed (Equation 6).

Environmental Conditions

Air temperature and air pressure on the three test days ranged from 12 to 13°C and 95.1–98.1 kPa, respectively. Average hourly wind speed from the two measurements stations ranged from 2.70 to 4.48 m·s⁻¹ (measured at 10 meters above ground), and average wind direction ranged from 14 to 46°.

Discussion

Accurate measurements of the propulsive power throughout a cross-country ski race could improve our understanding of the work requirements of cross-country skiing and other endurance sports that exhibit non-steady state power behavior (Figure 10). Calculations of propulsive power during cross-country skiing using the principle of power balance have been attempted by several authors (Moxnes and Hausken, 2008; Carlsson et al., 2011; Sandbakk et al., 2011; Moxnes et al., 2013, 2014; Swarén and Eriksson, 2017), but the accuracy of the method has not been thoroughly addressed. The current study addresses the accuracy obtained when applying the power balance principle to GNSS measurements during a roller skiing test race on a World Cup-like ski course. This was achieved using a Monte Carlo simulation sampling from the distributions of the most relevant sources of error (air drag, rolling resistance, variations in wind velocity), which enabled quantification of the method's accuracy throughout the course. Furthermore, the current study is the first to present continuous measurements of propulsive power throughout a test race in distance skiing (Sandbakk et al., 2011; Swarén and Eriksson, 2017). Our findings show that the error in the propulsive power estimates increases with skiing speed, while the propulsive power generated by the athletes decreases approximately linearly as speed increases. They also show that a substantial part of the variability of P_prop cannot be explained by skiing speed, accelerations, or measurement errors (Figure 5A). Explanations for this variability is most likely the complex course topography, which consisted of both relatively long uphill sections (i.e., from 2,000 to 2,500 m from the start, Figure 4D) and shorter uphill sections where the athletes had recovered during a long downhill section (3,600–3,730 m, Figure 4D). The duration of the longest uphill segment was typically ~120 s, while the short uphill from 3,600 to 3,730 m after start was completed in slightly more than 20 s. During all-out running or cycling tests with durations of 120 and 20 s, the anaerobic energy contributions are approximately 37 and 82%, respectively (Gastin, 2001). Therefore, it is a reasonable assumption that during the longer uphill sections P_prop is mainly limited by the athletes' VO_{2, max}, but this restriction does not apply to uphill sections with short durations, at least if the athletes are in a partially recovered state at the beginning of the uphill. This is in agreement with the observations in the current study, where P_prop appeared to converge to ~4 W·kg⁻¹ in the longest uphill while being almost 8 W·kg⁻¹ in the short uphill at 3,600–3,730 m after start.

FIGURE 10

Figure 10. (A) Heart rate (HR) and P_prop measurements during the third lap for an example athlete. HR is expressed relative to the peak HR measured during the test race. (B) Scatter plot showing all measurements of P_prop (normalized to each athlete's peak P_prop) vs. all measurements of HR. There was no significant correlation between the two parameters (r_Pearson = 0.01, p = 0.18). This indicates that changes in P_prop during cross-country ski races are too fast for heart rate to provide a valid measure of instantaneous metabolic power.

The observation that propulsive power is higher at low skiing speeds is in agreement with the notion that skiers focus their effort on the uphill sections, where the external resistance is increased due to a substantial component of gravity along the skiing surface. Uphill terrain is also known to be the terrain that is the major determinant of overall performance during time trials (Andersson et al., 2010; Sandbakk et al., 2011, 2016b; Bolger et al., 2015). This is consistent with conclusions in both cycling and cross-country skiing suggesting that athletes should increase their work rate in course segments where the external resistance is increased (Swain, 1997; Atkinson et al., 2007; Sundström et al., 2013). The rationale for this is that a decline in speed over a given distance is not compensated by an equivalent increase in speed over the same distance. This implies that athletes should to some extent aim at minimizing variations in speed by varying the propulsive power.

The peak power outputs measured in the current study are substantially below the values reported by Swarén and Eriksson (2017), who reported peak power outputs of 16 W·kg⁻¹ for a male skier during a classical style sprint race. Because the athlete analyzed by Swarén and Eriksson (2017) was a high-level skier (qualified for the final in a Continental Cup race), it would be fair to compare that skier to the best ranked skier in the current study. The best ranked skier in the current study had <15 FIS-points at the time of the data collection, and a peak power output of 8.6 W·kg⁻¹, which is still substantially below the findings of Swarén and Eriksson (2017). In contrast, our findings show a larger power output than the segment average power of 5.76 W·kg⁻¹ found by Sandbakk et al. (2011) during a sprint skating race. However, this was the average power output over relatively long uphill section (duration of ~1 min) and is therefore difficult to compare directly.

The relationship between propulsive power and metabolic power during ski skating is non-trivial (Sandbakk et al., 2012; Andersson et al., 2017); therefore we cannot deduce directly from our measurements how the metabolic energy demand depends upon skiing speed. However, if we use the measurement of gross efficiency at an 8% incline and a high work rate (η_Gross = 16.4 %) from Sandbakk et al. (2012) as a basis, our findings suggests a typical metabolic power of about 38 W·kg⁻¹ at skiing speeds of 2 m·s⁻¹, which corresponds to an oxygen demand of about 108 mL·kg⁻¹·min⁻¹ (Péronnet and Massicotte, 1991). For the peak power output measured, the corresponding oxygen demand would be 147 mL·kg⁻¹·min⁻¹. Assuming a $\dot{\text{V}}$O₂,_peak for this level of skiers of 75 mL·kg⁻¹·min⁻¹, this corresponds to a typical oxygen demand of about 144% of peak aerobic power at 2 m·s⁻¹, and a peak oxygen demand of about 196% of peak aerobic power, and must therefore elicit substantial use of anaerobic energy pathways. This is in agreement with previous studies in both sprint skiing (Sandbakk et al., 2011) and distance skiing (Norman and Komi, 1987; Karlsson et al., 2018). It is interesting to note that the peak metabolic energy requirements in sprint skiing and distance skiing appear to be relatively similar, even though the competition duration is substantially different (~2–4 min for sprint skiing and >30 min for a 15 km race). This might partly explain why cross-country skiing requires a relatively small degree of specialization to each discipline, allowing individual athletes to be world-class over distances ranging from ~1.3 to 50 km.

Methodological Considerations

The method's applicability to discriminate between P_prop and F_prop generated by high-level athletes depends on whether the measurement error is smaller than typical inter-athlete and intra-athlete differences observed throughout a race. Our findings show that the proposed method was not sufficiently accurate to discriminate between the inter-athlete or intra-athlete differences observed in this group of high-level skiers. The sources of these errors are distributed almost evenly between air drag (0.034 or 0.016 N·kg⁻¹ with measured wind and zero-wind, respectively) and rolling resistance (0.023 N·kg⁻¹). Because of varying surface-properties (i.e., asphalt quality) and the significant effect of temperature on the rolling resistance (Ainegren et al., 2008), it is challenging to obtain substantially more accurate measurements of rolling resistance. However, changes in rolling resistance caused by changes in asphalt quality or temperature will to some extent be systematic effects, and will therefore partially cancel when comparing intra- or inter-athlete differences during one experiment. Hence, the sensitivity-criterion used in the current study is appropriate when comparing results from different experiments, but might be too conservative for differences observed during a single experiment.

The results of this study clearly indicate the importance of precise measurements of vertical position. Estimates of propulsive power using vertical position measurements from the standalone GNSS receiver were substantially different (RMS deviation 31% of mean active propulsive power) from the measurements that were mapped onto the dGNSS reference trajectory. Hence, mapping standalone GNSS data on a precisely measured reference trajectory, or using athlete-mounted carrier phase differential GNSS receivers (Gilgien et al., 2014b; Karlsson et al., 2018) is required to calculate meaningful propulsive power measurements. Another solution that might be applicable is to fuse GNSS with IMU or barometric measurements (Skaloud and Limpach, 2003).

Predicting the drag area using only body mass and skiing speed yielded relatively small deviations from the accelerometer-based drag area model, and could therefore be an acceptable solution for many practical applications.

Limitations

In the current study we used roller skiing as an analog exercise for cross-country skiing on snow. Furthermore, a test race was used rather than an official ski race. Hence, the current study has lower external validity than studies performed on snow (Sandbakk et al., 2011; Swarén and Eriksson, 2017), but has higher internal validity because rolling resistance could be measured more accurately than ski/snow friction, and was less likely to change substantially during the test race.

Another challenge when applying the power balance principle is to accurately calculate air drag, because the drag area and wind speed are difficult to measure continuously along the course. As indicated in Figure 2, the drag coefficient of a cross-country skier depends on the Reynolds number. However, measurements are scarce and inconsistent, particularly at Reynolds numbers <2·10⁵. In the current study we created a model for how the drag coefficient changes with Reynolds number based on measurements on a brass cylinder (Achenbach, 1968), and scaled it to fit previous measurements of ice skaters at 12 m·s⁻¹ (van Ingen Schenau et al., 1982). We chose not to base our model on the measurements in Spring et al. (1988), as their equations did not account for changes in gravitational potential energy that would be caused by a non-level rolling surface. This could lead to substantial errors, particularly at low skiing speeds, even if the inclination of the rolling surface is very small. As an example, a 0.1° incline at 3 m·s⁻¹ would result in an error in C_DA of ~0.25 m². Hence, wind-tunnel studies investigating how the drag coefficient of a cross-country skier depends upon the Reynolds number would be useful, particularly at conditions relevant for low skiing speeds.

A challenge that was not addressed in the drag coefficient model proposed in current study is that cross-country skiing techniques causes body segments to move with different speeds through the air. This effect is most pronounced for the ski poles, where the pole tip's speed varies from 0 (when in contact with the ground) to an unknown speed substantially higher than the skiing speed. However, as the cross-sectional area of ski poles are relatively small, it is likely that this has only a minor effect on the total drag area.

Wind velocity was not measured continuously along the course but was estimated by a wind field based on the hourly average of two nearby meteorological stations, and the assumption of a Rayleigh distribution. The Monte Carlo approach simulated a large number (2,500) of different wind conditions and returned the expected value from all the simulated conditions, which should improve the calculations with respect to the assumption of a constant wind field. Nonetheless, it is obvious that the instantaneous wind velocity is strongly affected by gusts and the proximate surroundings of the course, and our calculations are therefore susceptible to such errors. The errors should however be within the uncertainty bounds specified by the Monte Carlo simulation.

The point mass assumption used in the current study neglects the work required to move body segments with respect to the athlete's center of mass, and therefore does not represent the total mechanical work done by the athlete. It is likely that the total mechanical work substantially exceeds the work required to move the center of mass alone. However, calculation of total mechanical work requires measurements of both the moments of force in all joints and the body segment kinematics (Aleshinsky, 1986; van Ingen Schenau and Cavanagh, 1990). Such measurements are currently not possible, at least in a field situation. Furthermore, even if the total mechanical work could be measured, there is no theoretically valid method to calculate the metabolic energy requirements based on kinetics and kinematics alone. Nevertheless, this should not discourage the use of statistical models linking propulsive power or total mechanical power to metabolic power under the assumption that the model is proven accurate for the problem being studied.

Two additional sources of error that have not been assessed in this study are that the skis do not move the exact same distance as the center of mass (due to the ski skating technique) and the fact that a small fraction of the surface normal force is exerted through the poles (estimated to 3–5% of body weight Millet et al., 1998). These effects have been assessed in other studies (Losnegard et al., 2012; Sandbakk et al., 2012) and are considered to be only of minor consequence.

Prospects

Studies in other winter sports where equipment-snow/ice friction and air drag are the main opposing forces have shown that the derivations of power, energy/work and propulsive force from athletes using position data are powerful approaches to studying the underlying mechanisms of performance (van Ingen Schenau et al., 1982; Supej, 2008; Gilgien et al., 2014a, 2016, 2018; Kröll et al., 2016). In endurance sports such as cross-country skiing, combining measurements of propulsive power with a model for skiing efficiency is a natural extension of the current study, and would improve our understanding of the physiological requirements of cross-country skiing (Sandbakk et al., 2011; Karlsson et al., 2018). Furthermore, simultaneous measurements of oxygen consumption would allow assessments of the balance between aerobic and anaerobic energy pathways at (or close to) the limit of human endurance racing performance (Andersson et al., 2017). Such measurements could also be used to improve numerical simulations of cross-country skiing performance (Moxnes et al., 2013, 2014), and to explore the effect of different pacing strategies (Swain, 1997; Sundström et al., 2013; Karlsson et al., 2018).

Although the Monte Carlo simulations used in the current study provides some insight into the method's validity and limitations, a validation against a gold standard has not been performed. Future studies could investigate the methods accuracy using ski poles and roller skis instrumented with force transducers.

Conclusion

During a 13.5 km roller skiing test race on a course similar to a cross-country skiing World Cup competition course, elite cross-country skiers generated a propulsive power output that declined approximately linearly with skiing speed, starting from 6.2 W·kg⁻¹ at the lowest measured speeds of 2.0 m·s⁻¹ (occurring at inclinations > ~ 10°). At skiing speeds close to 9 m·s⁻¹ and inclinations < ~ −2° the skiers transitioned to the tucked position where no propulsive power was generated. Furthermore, the results of this study clearly indicate the importance of precise measurements of vertical position, and shows that standalone GNSS receivers are not sufficiently accurate to be used for propulsive power calculations unless the measurements are mapped on a precisely measured reference trajectory, or replaced by carrier phase differential GNSS receivers. In contrast, predictions of drag area using only body mass and skiing speed deviated only slightly from those based on accelerometer data and should be acceptable for many practical applications. However, none of the methods presented in the current study were sufficiently accurate to discriminate between the instantaneous differences in propulsive force in this group of high-level athletes.

Author Contributions

Conception and design: ØG, TL, MG, DD, and AM-S. Data collection and data analysis: ØG and MG. Manuscript draft of Introduction: ØG and MG. Manuscript draft of Methods, Results, and Discussion: ØG. All authors contributed to manuscript revision, and read and approved the submitted version.

Funding

This study was funded by the Norwegian School of Sports Sciences, Oslo, Norway, and the Norwegian Research Council (project 216699).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

Achenbach, E. (1968). Distribution of local pressure and skin friction around a circular cylinder in cross-flow up to Re = 5^*10⁶. J. Fluid Mech 34, 625–639. doi: 10.1017/S0022112068002120

CrossRef Full Text | Google Scholar

Ainegren, M., Carlsson, P., and Tinnsten, M. (2008). Rolling resistance for treadmill roller skiing. Sport. Eng. 11, 23–29. doi: 10.1007/s12283-008-0004-1

CrossRef Full Text | Google Scholar

Aleshinsky, S. Y. (1986). An energy ‘sources' and ‘fractions' approach to the mechanical energy expenditure problem. J. Biomech. 19, 287–315.

Google Scholar

Andersson, E., Björklund, G., Holmberg, H.-C., and Ørtenblad, N. (2017). Energy system contributions and determinants of performance in sprint cross-country skiing. Scand. J. Med. Sci. Sports 27, 385–398. doi: 10.1111/sms.12666

PubMed Abstract | CrossRef Full Text | Google Scholar

Andersson, E., Supej, M., Sandbakk, Ø., Sperlich, B., Stoggl, T., and Holmberg, H. C. (2010). Analysis of sprint cross-country skiing using a differential global navigation satellite system. Eur. J. Appl. Physiol. 110, 585–595. doi: 10.1007/s00421-010-1535-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Atkinson, G., Peacock, O., and Passfield, L. (2007). Variable versus constant power strategies during cycling time-trials: prediction of time savings using an up-to-date mathematical model. J. Sports Sci. 25, 1001–1009. doi: 10.1080/02640410600944709

PubMed Abstract | CrossRef Full Text | Google Scholar

Bergh, U. (1987). The influence of body mass in cross-country skiing. Med. Sci. Sports Exerc. 19, 324–331.

PubMed Abstract | Google Scholar

Bolger, C. M., Kocbach, J., Hegge, A. M., and Sandbakk, Ø. (2015). Speed and heart-rate profiles in skating and classical cross-country skiing competitions. Int. J. Sports Physiol. Perform. 10, 873–880. doi: 10.1123/ijspp.2014-0335

PubMed Abstract | CrossRef Full Text | Google Scholar

Bowden, F. P., and Hughes, T. P. (1939). The mechanism of sliding on ice and snow. Proc. R. Soc. A Math. Phys. Eng. Sci. 172, 280–298. doi: 10.1098/rspa.1939.0104

CrossRef Full Text | Google Scholar

Buhl, D., Fauve, M., and Rhyner, H. (2001). The kinetic friction of polyethylen on snow: the influence of the snow temperature and the load. Cold Reg. Sci. Technol. 33, 133–140. doi: 10.1016/S0165-232X(01)00034-9

CrossRef Full Text | Google Scholar

Canuto, C., Quarteroni, A., Hussaini, M. Y., and Zang, T. A. (2007). Fundamentals of Fluid Dynamics. New York, NY: Springer.

Google Scholar

Carlsson, P., Tinnsten, M., and Ainegren, M. (2011). Numerical simulation of cross-country skiing. Comput. Methods Biomech. Biomed. Engin. 14, 741–746. doi: 10.1080/10255842.2010.493885

PubMed Abstract | CrossRef Full Text | Google Scholar

Dooner, D. B. (2012). Kinematic Geometry of Gearing. Chichester, UK: John Wiley & Sons

Google Scholar

Dorvlo, A. S. S. (2002). Estimating wind speed distribution. Energy Convers. Manag. 43, 2311–2318. doi: 10.1016/S0196-8904(01)00182-0

CrossRef Full Text | Google Scholar

Ettema, G., and Lorås, H. W. (2009). Efficiency in cycling: a review. Eur. J. Appl. Physiol. 106, 1–14. doi: 10.1007/s00421-009-1008-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Gastin, P. B. (2001). Energy system interaction and relative contribution during maximal exercise. Sport. Med. 31, 725–741. doi: 10.2165/00007256-200131100-00003

PubMed Abstract | CrossRef Full Text | Google Scholar

Gilat, A., and Subramaniam, V. (2008). Numerical Methods for Engineers and Scientists, 3rd Edn. Hoboken: NJ: John Wiley & Sons.

Google Scholar

Gilgien, M., Kröll, J., Spörri, J., Crivelli, P., and Müller, E. (2018). Application of dGNSS in alpine ski racing: Basis for evaluating physical demands and safety. Front. Physiol. 9:145. doi: 10.3389/fphys.2018.00145

PubMed Abstract | CrossRef Full Text | Google Scholar

Gilgien, M., Spörri, J., Kröll, J., Crivelli, P., and Müller, E. (2014a). Mechanics of turning and jumping and skier speed are associated with injury risk in men's World Cup alpine skiing: a comparison between the competition disciplines. Br. J. Sports Med. 48, 724–747. doi: 10.1136/bjsports-2013-092994

PubMed Abstract | CrossRef Full Text | Google Scholar

Gilgien, M., Spörri, J., Kröll, J., and Müller, E. (2016). Effect of ski geometry and standing height on kinetic energy: equipment designed to reduce risk of severe traumatic injuries in alpine downhill ski racing. Br. J. Sports Med. 50, 8–13. doi: 10.1136/bjsports-2015-095465

PubMed Abstract | CrossRef Full Text | Google Scholar

Gilgien, M., Spörri, J., Limpach, P., Geiger, A., and Müller, E. (2014b). The effect of different global navigation satellite system methods on positioning accuracy in elite alpine skiing. Sensors 14, 18433–18453. doi: 10.3390/s141018433

PubMed Abstract | CrossRef Full Text | Google Scholar

Günther, B. (1975). Dimensional analysis and theory of biological similarity. Physiol. Rev. 55, 659–699. doi: 10.1152/physrev.1975.55.4.659

PubMed Abstract | CrossRef Full Text | Google Scholar

Haugen, T., Paulsen, G., Seiler, S., and Sandbakk, Ø. (2017). New records in human power. Int. J. Sports Physiol. Perform. 13, 678–686. doi: 10.1123/ijspp.2017-0441

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoffman, M. D., Clifford, P. S., Bota, B., Mandli, M., and Jones, G. M. (1990). Influence of body mass on energy cost of roller skiing. Int. J. Sport Biomech. 6, 374–385. doi: 10.1123/ijsb.6.4.374

CrossRef Full Text | Google Scholar

Ingjer, F. (1991). Maximal oxygen uptake as a predictor of performance ability in women and men elite cross-country skiers. Scand. J. Med. Sci. Sports 1, 25–30. doi: 10.1111/j.1600-0838.1991.tb00267.x

CrossRef Full Text | Google Scholar

Karlsson, Ø., Gilgien, M., Gløersen, Ø., Rud, B., and Losnegard, T. (2018). Exercise intensity during cross-country skiing described by oxygen demands in flat and uphill terrain. Front. Physiol. 9:846. doi: 10.3389/fphys.2018.00846

PubMed Abstract | CrossRef Full Text | Google Scholar

Kröll, J., Spörri, J., Gilgien, M., Schwameder, H., and Müller, E. (2016). Sidecut radius and kinetic energy: equipment designed to reduce risk of severe traumatic knee injuries in alpine giant slalom ski racing. Br. J. Sports Med. 50, 26–31. doi: 10.1136/bjsports-2015-095463

PubMed Abstract | CrossRef Full Text | Google Scholar

Losnegard, T., and Hallén, J. (2014). Physiological differences between sprint- and distance-specialized cross-country skiers. Int. J. Sports Physiol. Perform. 9, 25–31. doi: 10.1123/ijspp.2013-0066

PubMed Abstract | CrossRef Full Text | Google Scholar

Losnegard, T., Myklebust, H., and Hallén, J. (2012). Anaerobic capacity as a determinant of performance in sprint skiing. Med. Sci. Sports Exerc. 44, 673–681. doi: 10.1249/MSS.0b013e3182388684

PubMed Abstract | CrossRef Full Text | Google Scholar

Millet, G. Y., Boissiere, D., and Candau, R. (2003). Energy cost of different skating techniques in cross-country skiing. J. Sports Sci. 21, 3–11. doi: 10.1080/0264041031000070903

PubMed Abstract | CrossRef Full Text | Google Scholar

Millet, G. Y., Hoffman, M. D., Candau, R. B., and Clifford, P. S. (1998). Poling forces during roller skiing: effects of technique and speed. Med. Sci. Sports Exerc. 30, 1645–1653. doi: 10.1097/00005768-199811000-00014

PubMed Abstract | CrossRef Full Text | Google Scholar

Moxnes, J. F., and Hausken, K. (2008). Cross-country skiing motion equations, locomotive forces and mass scaling laws. Math. Comput. Model. Dyn. Syst. 14, 535–569. doi: 10.1080/13873950802164643

CrossRef Full Text | Google Scholar

Moxnes, J. F., Sandbakk, Ø., and Hausken, K. (2013). A simulation of cross-country skiing on varying terrain by using a mathematical power balance model. Open Access J. Sport. Med. 4, 127–139. doi: 10.2147/OAJSM.S39843

PubMed Abstract | CrossRef Full Text | Google Scholar

Moxnes, J. F., Sandbakk, Ø., and Hausken, K. (2014). Using the power balance model to simulate cross-country skiing on varying terrain. Open access J. Sport. Med. 5, 89–98. doi: 10.2147/OAJSM.S53503

PubMed Abstract | CrossRef Full Text | Google Scholar

Myklebust, H. (2016). Quantification of Movement Patterns in Cross-Country Skiing Using Inertial Measurement Units. Dr. Thesis, Norwegian School of Sport Sciences.

Myklebust, H., Gløersen, Ø., and Hallén, J. (2015). Validity of ski skating center of mass displacement measured by a single inertial measurement unit. J. Appl. Biomech. 31, 492–498. doi: 10.1123/jab.2015-0081

PubMed Abstract | CrossRef Full Text | Google Scholar

Norman, R. W., and Komi, P. V. (1987). Mechanical energetics of world class cross-country skiing. Int. J. Sport Biomech. 3, 353–369. doi: 10.1123/ijsb.3.4.353

CrossRef Full Text | Google Scholar

Oke, T. R. (1978). Boundary Layer Climates, 2nd Edn. London: Methuen.

Google Scholar

Péronnet, F., and Massicotte, D. (1991). Table of nonprotein respiratory quotient : an update. Can. J. Sport Sci. 16, 23–29.

PubMed Abstract | Google Scholar

Sandbakk, Ø., Ettema, G., and Holmberg, H.-C. (2012). The influence of incline and speed on work rate, gross efficiency and kinematics of roller ski skating. Eur. J. Appl. Physiol. 112, 2829–2838. doi: 10.1007/s00421-011-2261-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Sandbakk, Ø., Ettema, G., Leirdal, S., Jakobsen, V., and Holmberg, H.-C. (2011). Analysis of a sprint ski race and associated laboratory determinants of world-class performance. Eur. J. Appl. Physiol. 111, 947–957. doi: 10.1007/s00421-010-1719-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Sandbakk, Ø., Hegge, A. M., Losnegard, T., Skattebo, Ø., Tønnessen, E., and Holmberg, H.-C. (2016a). The physiological capacity of the world's highest ranked female cross-country skiers. Med. Sci. Sports Exerc. 48, 1091–1100. doi: 10.1249/MSS.0000000000000862

PubMed Abstract | CrossRef Full Text | Google Scholar

Sandbakk, Ø., and Holmberg, H.-C. (2014). A Reappraisal of Success factors for olympic cross-country skiing. Int. J. Sports Physiol. Perform. 9, 117–121. doi: 10.1123/ijspp.2013-0373

PubMed Abstract | CrossRef Full Text | Google Scholar

Sandbakk, Ø., Losnegard, T., Skattebo, Ø., Hegge, A. M., Tønnessen, E., and Kocbach, J. (2016b). Analysis of classical time-trial performance and technique-specific physiological determinants in elite female cross-country skiers. Front. Physiol. 7:326. doi: 10.3389/fphys.2016.00326

PubMed Abstract | CrossRef Full Text | Google Scholar

Skaloud, J., and Limpach, P. (2003). “Synergy of CP-DGPS, accelerometry and magnetic sensors for precise trajectography in ski Racing,” in Proceedings of the 16th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GPS/GNSS) (Portland, OR).

Google Scholar

Spring, E., Savolainen, S., Erkkilä, J., Hämäläinen, T., and Pihkala, P. (1988). Drag area of a cross-country skier. Int. J. Sport Biomech. 4, 103–113. doi: 10.1123/ijsb.4.2.103

CrossRef Full Text | Google Scholar

Spurk, H., and Aksel, N. (2008). Fluid Mechanics, 2 Edn. Berlin: Springer-Verlag.

Google Scholar

Sundström, D., Carlsson, P., Ståhl, F., and Tinnsten, M. (2013). Numerical optimization of pacing strategy in cross-country skiing. Struct. Multidiscip. Optim. 47, 943–950. doi: 10.1007/s00158-012-0856-7

CrossRef Full Text | Google Scholar

Supej, M. (2008). Differential specific mechanical energy as a quality parameter in racing alpine skiing. J. Appl. Biomech. 24, 121–129. doi: 10.1123/jab.24.2.121

PubMed Abstract | CrossRef Full Text | Google Scholar

Swain, D. P. (1997). A model for optimizing cycling performance by varying power on hills and in wind. Med. Sci. Sports Exerc. 29, 1104–1108.

PubMed Abstract | Google Scholar

Swarén, M., and Eriksson, A. (2017). Power and pacing calculations based on real-time locating data from a cross-country skiing sprint race. Sport. Biomech. 3141, 1–12. doi: 10.1080/14763141.2017.1391323

CrossRef Full Text | Google Scholar

Theile, T., Szabo, D., Luthi, A., Rhyner, H., and Schneebeli, M. (2009). Mechanics of the ski–snow contact. Tribol. Lett. 36, 223–231. doi: 10.1007/s11249-009-9476-9

CrossRef Full Text | Google Scholar

van Ingen Schenau, G. J., and Cavanagh, P. R. (1990). Power equations in endurance sports. J. Biomech. 23, 865–881. doi: 10.1016/0021-9290(90)90352-4

PubMed Abstract | CrossRef Full Text | Google Scholar

van Ingen Schenau, G. J., Ohtsuki, T., Cruz, J. C., Hogan, R. P., and Balke, B. (1982). The influence of air friction in speed skating. J. Biomech. 15, 449–458. doi: 10.1016/0021-9290(82)90081-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Welde, B., Evertsen, F., Von Heimburg, E., and Ingulf Medbø, J. (2003). Energy cost of free technique and classical cross-country skiing at racing speeds. Med. Sci. Sport. Exerc. 35, 818–825. doi: 10.1249/01.MSS.0000064936.04725.FD

PubMed Abstract | CrossRef Full Text | Google Scholar