Frédéric Doumenc

Béatrice Guerrier

A model coupling the liquid and gas phases for a totally wetting evaporative meniscus

An hydrodynamic model has been developed to get a complete description of an evaporative meniscus in complete wetting configuration. The coupling between the liquid and gas is explicitly taken into account. Scaling laws are derived for the different domains of the meniscus and validated by numerical simulations. Results are compared with previous models of the literature that use the electrostatic analogy proposed by Deegan and co-authors to describe the evaporative flux. We show that the different approaches differ for the description of the tip of the meniscus in the domain corresponding to the decrease of the evaporative flux but lead to the same scaling for the apparent macroscopic contact angle.

Introduction

Evaporation of droplets or meniscus has given rise to many studies in the last years, both in terms of fundamental understanding [1] and of processing of patterned deposits (for instance [2][3][4]). Local description of the contact line for an evaporative droplet or meniscus is still an open question and several approaches have been recently developed to take into account the coupling between evaporation and hydrodynamics. An important distinction to be drawn concerns the gas phase surrounding the evaporating liquid. Indeed the phenomena involved are very different depending on wether the liquid is surrounded by its own vapor only [5][6][7] or by an inert gas. In the first configuration evaporation is driven by the departure from equilibrium at the interface while, when evaporating in an inert gas, the limiting phenomenon comes from the diffusion of the vapor in the gas [8]. In this paper we restrict the discussion to the second configuration, i.e. evaporation in an inert gas.

More specifically we focus on total wetting configurations where evaporation can modify the behavior observed for a non volatile fluid. For instance, Cazabat and coauthors [9][10][11] have performed numerous experiments on volatile droplets left on a substrate and observed the spreading of the droplet up to a maximum radius, followed by retraction due to the competition between wetting and evaporation. For meniscus on a moving substrate (like dip coating experiments), evaporation also changes significantly the meniscus shape when the substrate velocity is not too large (for large substrate velocities, classical Landau-Levich or Cox-Voinov models apply).

a e-mail: doumenc@fast.u-psud.fr

Finally we restrict the study to pure fluids. Indeed, as far as complex fluids are concerned, other hydrodynamic behavior of the tip is expected due to the great change of the viscosity with solute concentration. This has been investigated in a previous study where drying of polymer solution was analyzed [12].

In this framework (evaporation in an inert gas, total wetting, pure fluid) this paper presents a numerical study and scaling laws for the different regions of the meniscus. As noticed by Deegan and co-authors [13,14], an electrostatic analogy allows to derive an estimation of the profile of the evaporative flux, as far as local description in the vicinity of the drop edge is not required. Indeed, using the electrostatic analogy, the evaporative flux is shown to diverge close to the droplet contact line. This divergence is of course not physical and one of the objective of this paper is to evaluate the validity of that model by comparing it with a more sophisticated approach including the coupling between the gas phase and the liquid. We thus compare our results to different models developed in the literature [9,15,16].

Problem formulation

Model equations

The geometry considered here corresponds to the capillary rise of a volatile fluid between two parallel plates immersed in an infinite reservoir. We consider a stationary configuration where the meniscus is motionless. The temperature is assumed uniform (isothermal problem). Coupled mass transfers in the liquid and in the gas phases are taken into account. Drying takes place in air at atmospheric pressure, local thermodynamic equilibrium is assumed at the interface and the evaporation is limited by the solvent vapor transfer into stagnant air (1.5-sided model).

In the liquid phase the velocity field is assumed quasi-parallel to the substrate (lubrication approximation) with no-slip at the substrate and zero tangential stress at the free surface. The direction of the flow is noted x, the axis normal to the substrate is z. Gravity is neglected.

This model is very close to the one developed previously in the same geometry for a binary solution on a moving substrate [12], or to the one used by Eggers and Pismen for a drop of volatile liquid [15].

With the above-mentioned approximations Stokes equations lead to the standard result [5]:

Q(x, t) = h 3 3η ∂ ∂x (σ ∂ 2 h ∂x 2 + Π d ) ( 1
)
where h is the liquid height, η the dynamic viscosity, σ the surface tension and Π d the disjoining pressure. To model total wetting we use a long-range stabilizing Van der Walls disjoining pressure Π d = A/(6πh 3 ) with the Hamaker constant A > 0. Q(x, t) is the liquid volumic flux (by unit of width) at abscissa x and time t over a cross section normal to x,

Q(x, t) = ∫ h(x) 0 u(x, z)dz (2)
where u(x, z) is the velocity component in the x direction. The mass conservation reads:

∂h ∂t + ∂Q ∂x = -v ev (x, t) (3)
where v ev (x, t) is the local evaporation velocity (positive for evaporation).

The boundary conditions are the following: at x = 0 we impose the liquid height h 0 and the curvature C 0 . If one performs an analogy with previous capillary rise experiments [18], C 0 is the curvature of a static meniscus, deduced from the gap d between the plates: C 0 = 2/d, with d = 1mm. The lubrication theory is valid if the free surface slope is much lower than 1, which is true if h 0 is small enough.

With the geometric assumptions of the model, we define a length L st used in the following as the order of magnitude of the meniscus size: considering a parabolic meniscus (constant curvature C 0 and small slope approximation) of height h 0 at x = 0, it intercepts the substrate at x = L st with zero contact angle. A simple derivation gives L st = √ 2h 0 /C 0 . Spatial derivatives of h in equation 1 being of fourth order, two boundary conditions are also needed at x = W . These conditions far from the bulk are given in equation 5 (as far as the length of the simulation box is large enough, we have checked that the choice of the outlet conditions does not affect the results).

at x = 0 : h = h 0 ∂ 2 h ∂x 2 = C 0 (4) at x = W : ∂h ∂x = 0 ∂ 2 h ∂x 2 = 0 (5)
The local evaporation velocity in equation 3 is determined by solving a diffusion problem in the gas phase, in a domain of length W and height H (for process with a forced air flow above the free surface, H would roughly correspond to the boundary layer thickness). Since H is much larger than h 0 , we assume a rectangular domain for the gas phase. This means that liquid thickness variations are neglected when computing mass transfer in the gas phase. The bottom boundary of this domain corresponding to the liquid-gas interface is then: z = 0 and 0 < x < W . The diffusion equation in the gas phase reads:

∂c g ∂t = D g ( ∂ 2 c g ∂x 2 + ∂ 2 c g ∂z 2 )
for 0 < x < W and 0 < z < H (6) where c g is the solvent vapor concentration in the air. The vertical walls are assumed impermeable. At the top a Dirichlet condition is imposed, with solvent vapor concentration equal to c gH . Then we get the following boundary conditions for the gas domain:

∂c g ∂x = 0
for (x = 0 or x = W ) and 0 < z < H (7)

c g = c gH for 0 ≤ x ≤ W and z = H (8)
The last boundary condition corresponds to the coupling between the liquid and the gas phases and is achieved by writing the mass flux conservation and the local thermodynamic equilibrium at the interface:

ρ v ev = -D g ∂c g ∂z for 0 ≤ x ≤ W and z = 0 ( 9
)
c g = c gs0 exp[- vs RT (σ ∂ 2 h ∂x 2 + Π d )] for 0 ≤ x ≤ W and z = 0 ( 10
)
with c gs0 = Ms RT P vs0 (T ). ρ is the liquid solvent density, M s the liquid molar mass, R the ideal gas constant, T the temperature, P vso the bulk saturated vapor pressure of the fluid, and vs its molar volume. As already said, evaporation cooling effect is disregarded so the temperature T is assumed uniform in the whole system, and equal to the surrounding temperature.

This paper focuses on the steady state and we use as initial conditions the result of a previous calculus for a close configuration. It has been checked that the choice of the initial conditions has no effect on the steady state result.

The above set of equations (1 to 10) is solved numerically by finite differences, using a pure implicit scheme of order 1 in time and 2 in space. At each time step equations are solved iteratively in the liquid and in the gas, until convergence is achieved. The discretization of equations in the liquid (1 to 5 and 9) leads to a set of nonlinear algebraic equations, solved by Newton-Raphson algorithm. An adaptive mesh is implemented, in order to refine the mesh in the region of high evaporative flux. In the gas phase, the set of algebraic equations resulting from discretization of equations 6 to 8 and 10 is linear, and solved by Successive Over-Relaxation method (SOR). A specific software has been written in Fortran, and most of the computations have been performed on iDataPlex IBM clusters. Computational times vary from a few hours to several weeks for high evaporative flux.

Results

Local description

We first present a detailed analysis of the meniscus shape, the liquid flow, the local evaporative flux and the concentration field in the gas for a reference configuration corresponding to the following data: the box dimensions are W = 1mm and H = 3mm, the boundary conditions are h 0 = 0.1mm, C 0 = 2mm -1 and c gH = 0. Physical properties are close to properties of toluene at 25

0 C: ρ = 900kg.m -3 , σ = 30mN.m -1 , M s = 92.14g.mol -1 , η = 0.55mP a.s.
Vapor diffusivity in the air is D g = 8.6 × 10 -6 m 2 .s -1 and the saturated vapor pressure is given by Antoine equation: log(P vs0 ) = a 1 -a 2 /(T + a 3 ), with a 1 = 9.0782, a 2 = 1343.9, a 3 = -53.77, P vs0 in Pascal and T in Kelvin. The Hamaker constant is A = 10 -19 J.

The profile of the meniscus is given in figure 1(a). From the observation of the meniscus slope (figure 1(b)), the liquid pressure (figure 1(c)) and the evaporation rate (figure 1(e)), we can split the meniscus into four characteristic domains, with the corresponding boundaries x 1 to x 4 :

-Bulk region (0 ≤ x ≤ x 1 ). In the bulk the slope ∂h/∂x varies linearly with x (i.e.

quasi constant curvature approximately equal to C 0 ) and the meniscus profile is then parabolic. The first boundary x 1 corresponds to a significant change in the meniscus curvature, arbitrarily set to 10%. In the following we call macroscopic angle the value of the slope at x = x 1 . -Capillary domain (x 1 ≤ x ≤ x 2 ). This domain is characterized by a strong increase of the curvature. The pressure in the liquid film is dominated by the capillary term. The boundary x 2 corresponds to equality between the capillary pressure and disjoining pressure. -Precursor film (x 2 ≤ x ≤ x 4 ). In this domain the disjoining pressure dominates over the capillary pressure and the evaporative flux is non negligible. It can be split into two sub-domains: a region of increasing evaporative flux (x 2 ≤ x ≤ x 3 ) followed by a region of decreasing evaporative flux (x 3 ≤ x ≤ x 4 ). The evaporative flux is thus maximum for x = x 3 . -Adsorbed film (x > x 4 ). Finally the meniscus ends by a film of nanometric thickness where the evaporative flux and liquid flow are almost zero. We call x 4 the abscissa where the evaporative flux is less than 3% of its maximum value. Evaporation is assumed negligible for x > x 4 .

The four boundaries characterizing the different domains of the meniscus (based on order of magnitude arguments) are reported on all the figures describing the local behavior of the meniscus. The local shape of the meniscus results from the coupling with the 2D vapor concentration field in the gas phase. Indeed a strong deformation of isoconcentration lines in the gas phase is observed at the meniscus tip, as shown in figure 1(d).

Evaporative flux

In the following we often argue about orders of magnitude. In that case the numerical prefactors are disregarded. The total evaporative flux I and the mean evaporation rate vev read:

I = ∫ W 0 v ev (x)dx and vev = 1 x 4 ∫ W 0 v ev (x)dx ≃ 1 L st ∫ W 0 v ev (x)dx (11)
For x 4 = W (length of the box), the vapor transfer in the gas is 1D and the total evaporative flux simply reads

I 1D = D g ρH (c gs0 -c gH )W
For x 4 < W and W << H, the isovalues of the vapor concentration are similar to those of the 1D configuration in most part of the box and the total evaporative flux is then almost the same (cf [12] for more details). Then we get the following scalings

vev ∼ v sc ev with v sc ev ≡ D g ρH (c gs0 -c gH ) W L st ( 12
)
This scaling is still relevant in our case where H = 3W , as shown in figure B.1(a). For a droplet in a semi-infinite medium, one should choose H ∼ W ∼ L st ∼ R, with R the droplet radius.

As already noticed in the introduction, Deegan and co-authors use an electrostatic analogy to get an estimation of the profile of the evaporative flux for a droplet, which scales as 1/ √ x div -x in the vicinity of the droplet's edge and diverges at the droplet's edge x div . Several authors use this expression to model the evaporation of droplets or to compute macroscopic contact angle. It is then interesting to compare the evaporative flux deduced from numerical simulation to this law. In this aim of view, the evaporative rate deduced from the above model is approximated by w ev (x) = J 0 / √ x div -x. The two parameters J 0 and x div are fitted to get w ev (x) ≃ v ev (x) for x x 3 (i.e. before the peak of the evaporative flux). For the reference configuration we get J 0 = 9.78.10 -9 m 3/2 s -1 and x div = 293.56µm. The two evaporative fluxes v ev (x) (numerical simulation) and w ev (x) (electrostatic analogy) are compared in figure 1(e) for the reference case. As can be seen the approximation is consistent with our model for x x 3 . At the very end of the meniscus evaporation fades away and then the electrostatic analogy does not hold anymore. Figure 1(f) shows the integrals of the evaporative flux from the end of the box (x = W ) for the two expressions v ev (x) and w ev (x), that is: Since the box is large enough to get v ev (W ) ≃ 0, since the problem is stationary and the substrate is motionless, mass conservation implies the equality of I(x) with the hydrodynamic flux Q(x). On the other hand it can be seen that for x x 3 , i.e. before the peak of the evaporative flux, the two integrals I(x) and I w (x) are very close (figure 1(f)), which means that the hydrodynamic flux that goes through the film for x x 3 is the same for the two representations of the evaporative flux.

I(x) = ∫ W x v ev (x)
We are now going to derive the scaling laws for the maximum of the evaporative flux v max ev and the length (x 4 -x 3 ) corresponding to the decrease of the evaporative flux. For x ∼ x 3 an estimation of the integrals is derived from the two models.

The flux is maximum for x = x 3 and varies regularly until x ≃ x 4 , then we can write

Q(x 3 ) = I(x 3 ) ∼ (x 4 -x 3 ) v max ev with v max ev = v ev (x 3 ) (13)
I w (x 3 ) can be deduced directly from the analytical expression of w ev (x):

I w (x 3 ) ∼ J 0 √ x div -x 3
The numerical simulations suggest that I(x 3 ) ∼ I w (x 3 ) and v ev (x 3 ) ∼ w ev (x 3 ). Then we easily deduce the following relation:

(x 4 -x 3 ) ∼ (x div -x 3 ) (14)
At x = x 3 the capillary pressure is much lower than the disjoining pressure and can be disregarded so that equation 1 reads:

Q(x 3 ) ∼ A ηh ∂h ∂x (15)
In this domain the film is thin enough for the disjoining pressure to affect the saturated vapor pressure so that the evaporative flux and film height are strongly coupled.

We then assume that the evaporative flux and film height decrease of one order of magnitude on the same length (x 4 -x 3 ). Using I(0) ∼ I w (0) and equations 13, 14 and 15 lead to the following scaling:

(x 4 -x 3 ) ∼ (x div -x 3 ) ∼ ξ with ξ ≡ ( A η v sc ev √ L st ) 2/3 (16)
At last, approximating the total length of the meniscus x 4 by L st , we easily get:

v max ev ∼ v sc ev √ L st /ξ ( 17
)
Equations 16 and 17 are the desired scalings. It is worthwhile to notice that the two lengths corresponding to the tip of the meniscus, the real one (x 4 -x 3 ) and the fictive one (x div -x 3 ), both scale as ξ.

Scaling laws of local slopes and thicknesses

This section is dedicated to the derivation of the scaling laws allowing the local description of the meniscus for the different domains described in section 3.1. Moving back towards the bulk we express in the vicinity of x = x 2 the balance between the effect of the capillary pressure, the disjoining pressure and the evaporation. Using the

(x2 -x1) ∼ a 1/3 /C 2/3 0 h1 ∼ a 5/6 /(ξ 1/2 C 2/3 0 ) θ1 ∼ (a/ξ) 1/2 (x4 -x2) ∼ ξ h2 ∼ (aξ) 1/2 θ2 ∼ (a/ξ) 1/2 (x4 -x3) ∼ ξ h3 ∼ λ 2/3 ξ 1/3 θ3 ∼ (λ/ξ) 2/3
same arguments than in the previous section, the balance between the flux induced by the disjoining pressure and the evaporation leads to the same scaling, i.e.

(x 4 -x 2 ) ∼ ξ At x = x 2 the pressures are equal so that we get

A h 3 ∼ σ ∂ 2 h ∂x 2
Assuming that the height and slope decrease of at least one order of magnitude between x 2 and x 4 (cf figure 1), we obtain

h 4 2 ∼ ξ 2 A/σ and θ 2 ∼ h 2 /ξ that is h 2 ∼ √ ξa and θ 2 ∼ √ a/ξ ( 18
)
with a ≡ √ A/σ, h 2 ≡ h(x 2 ) and θ 2 the local slope at x = x 2 . Using the same approach, scaling laws have been derived for the film heights, h i ≡ h(x i ), the domain lengths and the local slopes in the different regions of the meniscus. They are summarized in Table 1 and detailed in the appendix A. For h 1 and θ 1 they are derived from the approximated expression of the flux, that is w ev (x). In the region such that x 3

x x 4 , the local description of the coupling between the liquid and the gas must be taken into account to get h 3 and θ 3 .

This analysis shows that the relevant length scales involved in the problem are the following:

λ ≡ vs RT √ ηP vs0 D g , a ≡ √ A σ , ξ ≡ ( A ηv sc ev √ L st ) 2/3 (19)
λ and a are microscopic lengths which depend only on the fluid physical properties. For most of the common volatile fluids at room temperature, λ ∼ 0.1nm and a ∼ 1nm.

Only the characteristic length ξ depends on the evaporative flux. In that study, it has been varied from 0.1 to 1000µm. The geometry of the problem appears through the length L st = √ 2h 0 /C 0 . In order to confirm the above scalings and their validity domain, simulations were performed for a large range of parameters, i.e. the humidity, temperature and system properties. The parameters of the simulations are resumed in Table 2. Results of numerical simulations compare very well with scaling laws, as shown in figures 2 and B.1. The limit of validity of the scalings can clearly be observed on the figures and are detailed in next section.

Validity domain of the scaling laws

A first limitation of the above description is obtained when the peak of the evaporative flux occurs in the domain where the capillary pressure is dominant, i.e. when x 3 becomes smaller than x 2 . Indeed the arguments used to derive ξ neglect the capillary pressure which is no more valid in this configuration. The first limit of validity is then given by the condition h 2 h 3 , that is

ξ 1/6 a 1/2 λ -2/3 1
However this must be viewed as a theoretical limit. Indeed, in realistic configurations of evaporation in an inert gas, it seems that no real liquid can go beyond this limit (even if the prefactors are taken into account).

The other limitation corresponds to the violation of the condition h 2 h 1 , which leads to: ξa -1/3 C 2/3 0 1 Beyond this limit evaporation turns to be negligible and the meniscus profile tends to the classical solution corresponding to thermodynamic equilibrium for a non volatile fluid. From pressure equality between the film and the bulk domain, we get

h 1 ∼ h 2 ∼ a 2/3 /C 1/3 0 .

Comparison with other models

In this last section our results are compared to the models proposed by Cazabat and co-authors [9][10][11] on one hand, and Pham and co-authors [16] on the other hand. The first model is dedicated to droplet evaporation, while the later addresses the problem of an evaporating wedge, with zero curvature far from the contact line. Total wetting is assumed in both cases, and the evaporative flux is described by Deegan's electrostatic analogy. The main difference in the two approaches lies in the treatment of the precursor film.

-Cazabat et al [10] assume that close to the droplet edge, the leading terms in equations 1 and 3 are the capillary and the disjoining pressures, and then neglect advection and evaporation in that region. This leads to a tip profile which scales like √ ax. The crossover to the macroscopic region with contact angle θ is found to take place at a distance of order aθ -2 , where the disjoining pressure turns to be negligible, and the evaporation rate reaches its maximum value. -Pham et al ( [16] and erratum [17]) fix a priori the height of the film and a zero slope at the tip of the meniscus (corresponding to our x = x 4 ) where they impose a height which is equal to our expression of h 2 (while our numerical simulations give a very different thickness: h(x 4 ) ∼ 1nm). The problem is solved by an approximated analytical method, validated by numerical simulations.

For both approaches, the estimation of the length of the precursor film (where disjoining pressure in non negligible) is found to be of order ξ, which is consistent with our result (x 4 -x 2 ) ∼ ξ. Moreover, in the limit of zero contact line velocity (capillary number Ca = 0), these two models give the same scaling for the apparent contact angle:

θ ∼ (η vev ) 1/3 L 1/6 A 1/12 σ 1/4 ( 20
)
where θ is the contact angle, vev is the mean evaporation rate, L is the droplet radius or the wedge size. It happens that this scaling is the same than the one we obtained for θ 1 or θ 2 (see Table 1). In addition, Pham et al find a prefactor equal to 2.5, which is very close to our value 2.7 for θ 1 (see figure 2c).

It is thus worthwhile to note that the three approaches lead to the same scaling for the macroscopic angle, while the description of the tip of the droplet or meniscus (for x > x 3 ) are very different. We suggest the following explanations:

the squareroot profile of Deegan's electrostatic analogy succeeds in giving the right integral of the evaporative flux outside the region where it falls down (i.e. for x x 3 ). -In all these models, there exists a domain of the droplet or the meniscus where fluxes corresponding to evaporation, capillary forces and Van der Waals forces are all of same order of magnitude. According to scaling analysis of section 3.3, this seems to be a key ingredient.

However models based on electrostatic analogy are not able to describe the end of the precursor film and it is necessary to shift to the complete model presented in this paper to get the film profile when the flux begins to decrease (x > x 3 ). Finally, Eggers and Pismen [15] have studied the dynamic of evaporative drops retraction and have developed a model similar to the one presented in this study, taking into account explicitly the coupling between the liquid and the gas. They got the same scaling ξ for the width of the contact line region.

Conclusion

In this paper an hydrodynamic model has been developed to get a complete description of an evaporative meniscus in the configuration of complete wetting. The coupling between the liquid and gas is explicitly taken into account. Scaling laws based on heuristic arguments are derived for the different domains of the meniscus and validated by numerical simulations. Results are compared with previous models of the literature using the electrostatic analogy suggested by Deegan and co-authors. We show that the different approaches differ for the description of the tip of the meniscus in the domain corresponding to the decrease of the evaporative flux but lead to the same scalings for the apparent macroscopic contact angle. A possible extension should be to release the assumption of isothermal state [19]. Finally this model is based on a continuum description of the fluid. For the description of the tip of the meniscus, at the scale of the junction of the meniscus to the adsorbed film of nanometric thickness, the pertinence of such a model may be questionable. However this should have no effect on the estimation of the macroscopic apparent contact angle.

A Scaling analysis A.1 Maximum of evaporation velocity

We already established in section 3.2 that the evaporation flux decreases over a length:

(x 4 -x 3 ) ∼ ξ with ξ ≡ ( A η v sc ev √ L st ) 2/3 (A.1)
We derive now the height h 3 and the angle θ 3 at x = x 3 , where the evaporation flux reaches its maximum. In the region of deacreasing evaporation flux, we must take into account the coupling between the liquid and the gas phases due to the effect of disjoining pressure on saturated vapor pressure (capillary pressure can be neglected in that region). After linearizing equation 10, we get the following scaling for the thickness:

h 3 ∼ ( vs A RT c gs0 ∆c g ) 1/3 (A.2)
with ∆c g the variation of vapor concentration in the gas phase above the liquid film, in the region of decreasing evaporation flux of characteristic length ξ. Assuming the same scaling for the horizontal and vertical components of concentration gradient (∂c g /∂x ∼ ∂c g /∂z), the boundary condition 9 gives a relation between ∆c g and v max ev :

D g ∆c g ξ ∼ ρv max ev (A.3)
Using equations A.2,A.3 along with 17 for v max ev , we get the desired scaling:

h 3 ∼ λ 2/3 ξ 1/3 with λ ≡ vs RT √ ηP vs0 D g (A.4)
The angle θ 3 is then obtained by writing θ 3 ∼ h 3 /ξ.

Fig. 1 .

 1 Fig. 1. Meniscus structure, (a)Liquid height (semilogarithmic scale, linear scale in the insert), (b)Slope, (c)Liquid pressure, (d)Normalized concentration isovalues in the gas phase cg/cgs0, (e)Evaporation rate, (f)Integral of evaporation rate.

Fig. 2 .

 2 Fig. 2. Validation of scaling laws. Symbols : results from numerical simulations (see Table 2), lines : fit (see legend). (a) Maximum of velocity, (b) Horizontal lengths, (c) Local slopes at x = x1 and x = x2, (d) Local slope at x = x3.

Fig. B. 1 .

 1 Fig. B.1. Validation of scaling laws. Symbols : results from numerical simulations (see Table 2), lines : fit (see legend). (a) Mean evaporation rate, (b) Film thickness at x = x1, (c) Film thickness at x = x2, (d) Film thickness at x = x3.

Table 1 .

 1 Scaling laws.

Table 2 .

 2 Parameters used in the simulations. Only the parameters different from the reference case are noted. Symbols are used in figures 2 and B.1. (*): Pvs0 computed from Antoine's law.

	Serie	Symbol	T	cgH /cgs0	Lst	Dg	Pvs0	A	η	σ
			0 C		µm	mm 2 .s -1	P a	J	mP a.s	mN.m -1
	Ref. case		25	0	316	8.6	(*)	10 -19	0.55	30
	#1	△		[0; 0.99]
	#2	∇		[0; 0.99]				10 -20
	#3	▹		0.5		[3; 100]
	#4	◃	[-10; 52]
	#5	♢					[1; 3792]
	#6	+		0.9					[0.1; 550]
	#7	×		0.5					[0.1; 55]
	#8	⋆	6.85						[0.1; 55]
	#9	•		0.5						[10, 100]
	#10	*	6.85							[10, 100]
	#11			0.5	[100; 316]

Acknowledgments

The authors gratefully acknowledge the financial support from the MULTIFLOW Marie Curie FP7-ITN Network and have benefited from fruitful discussions with H.Bodiguel, F.Lequeux, L.Limat and C. Pham. The numerical calculations have been performed on the computers of the "Direction des Services d'Information, Service Calcul" of the university UPMC.

A.2 Boundary between the bulk and the capillary zone

We see from numerical simulations that the slopes at x = x 1 and x = x 2 have the same scaling: θ 1 ∼ θ 2 . Let us denote l c ≡ (x 2 -x 1 ) the length of the capillary zone (liquid pressure dominated by the capillary term along with significant curvature gradient). The bulk region is characterized by a negligible curvature gradient, i.e. ∂C/∂x ≪ C 0 /l c . Therefore, at the boundary between bulk and capillary zones (x = x 1 ), we get the relation :

Keeping the leading terms in equations 1 and 3, and assuming l c ξ, we get:

Using equations A.5, A.6, 18 and relation I(0) ∼ I w (0) to get J 0 , we obtain the following scalings:

Violation of condition l c ξ (i.e. l c ≪ ξ) corresponds to the transition to another regime, in which evaporation is negligible (see section 3.4).

B Numerical simulations : complementary results

We present in figure B.1 some results of numerical simulations to complete the validation of scaling laws presented in Table 1


Ph.D Najoua Assila
email: najoua.assila@umontreal.ca

Ph.D Paula W Rushton

Sonia Duprey

Ph.D.; Mickaël Begon

Trunk and glenohumeral joint adaptations to manual wheelchair propulsion over a cross-slope: an exploratory

Keywords: Wheelchair propulsion, Cross-slope, Glenohumeral joint, Trunk, Wheelchair skills

Background: Cross-slopes are often encountered by manual wheelchair users propelling within an urban setting. While propulsion over cross-slopes is more difficult than on level surfaces, little is known about how the users counter the downhill turning tendency of the wheelchair over cross-slopes. This study aimed to identify the adaptations of the manual wheelchair users to the presence of cross-slopes and examine how these might impact shoulder injury.

Methods: Nine manual wheelchair users propelled themselves across a cross-slope and over a level surface. The trunk and glenohumeral joint kinematics, as well as the handrim contact tangential force were compared between both conditions for the uphill and downhill limbs.

Findings:

The uphill arm technique used to counter the downhill turning tendency varied greatly in terms of potential injury risk and efficiency between participants. Trunk flexion increased the turning tendency of the manual wheelchair, yet only one participant decreased his flexion when rolling over the cross-slope. Various potential pathomecanisms related to the trunk lateral flexion and the glenohumeral kinematics over a cross-slope were identified.

Interpretation: Both the uphill arm technique and trunk kinematics are important to propel over a cross-slope both efficiently and safely. Accordingly, tips about posture and kinematics are needed to teach this skill to manual wheelchair users. Additionally, as wheelchair positioning seems to influence the cross-slope skill, more research is needed to explore the impact of positioning devices (e.g., lateral supports) and wheelchair modifications (e.g., power assist wheels, handrim projections) on this skill.

INTRODUCTION

The quality of life and the health of manual wheelchair users (MWUs) are associated with their participation. [START_REF] Levasseur | Subjective Quality-of-Life Predictors for Older Adults with Physical Disabilities[END_REF] Pain and wheelchair skills are two frequently reported modifiable factors that influence social participation. [START_REF] Smith | A review of factors influencing participation in social and community activities for wheelchair users[END_REF][START_REF] Jensen | Symptom Burden in Persons With Spinal Cord Injury[END_REF][START_REF] Mortenson | Association Between Mobility, Participation, and Wheelchair-Related Factors in Long-Term Care Residents Who Use Wheelchairs as Their Primary Means of Mobility[END_REF] Within an urban setting, MWUs will frequently encounter cross-slopes. Indeed, roads and sidewalks alike are slightly sloped to help water drainage. As MWUs propel themselves over a cross-slope, their wheelchair will tend to turn downhill due to the moment generated by the weight of the manual wheelchair and its user. [START_REF] Brubaker | Effects of side slope on wheelchair performance[END_REF] While there is a consensus that cross-slope propulsion is more difficult than propulsion over a level surface, [START_REF] Cooper | Manual Wheelchair Propulsion Over Cross-Sloped Surfaces: A Literature Review[END_REF] little information on its impact on the MWU can be found in the literature. [START_REF] Cooper | Manual Wheelchair Propulsion Over Cross-Sloped Surfaces: A Literature Review[END_REF][START_REF] Rouvier | Manual wheelchair biomechanics while overcoming various environmental barriers: A systematic review[END_REF] Brubaker et al. found that the presence of a 2° cross-slope roughly doubled the drag at similar velocities, increasing the propulsion energy cost and oxygen consumption. [START_REF] Brubaker | Effects of side slope on wheelchair performance[END_REF] The propulsion speed seems to decrease, particularly due to the decrease in the distance traveled per push. [START_REF] Richter | Consequences of a Cross Slope on Wheelchair Handrim Biomechanics[END_REF] Peak handrim force, rate of loading, and axial moment measured for the downhill wheel were reported to increase with the cross-slope angle. [START_REF] Richter | Consequences of a Cross Slope on Wheelchair Handrim Biomechanics[END_REF] On the other hand, the uphill wheel was reported to have a negative mechanical power as it regulates the wheelchair trajectory. [START_REF] Sauret | Handrim mechanical power during wheelchair propulsion on level and cross-slope surfaces: a preliminary study[END_REF] As the amplitude of handrim forces, the rate of loading and shoulder kinematics are correlated with musculoskeletal injury, [START_REF] Mercer | Shoulder joint kinetics and pathology in manual wheelchair users[END_REF][START_REF] Briley | The longitudinal relationship between shoulder pain and altered wheelchair propulsion biomechanics of manual wheelchair users[END_REF][START_REF] Boninger | Shoulder magnetic resonance imaging abnormalities, wheelchair propulsion, and gender[END_REF] the increased asymmetry and biomechanical load [START_REF] Sauret | Handrim mechanical power during wheelchair propulsion on level and cross-slope surfaces: a preliminary study[END_REF][START_REF] Chénier | Proposing a new index to quantify instantaneous symmetry during manual wheelchair propulsion[END_REF] observed for cross-slope propulsion might increase musculoskeletal injury. Thus, this skill requires specific training of the MWUs to insure a safe and efficient technique. [START_REF] Kirby | Wheelchair Skills Program Manual Version 5.2[END_REF] Nevertheless, we could not find studies that gave recommendations about cross-slope propulsion technique.

The Wheelchair Skills Program [START_REF] Kirby | Wheelchair Skills Program Manual Version 5.2[END_REF] (WSP) includes cross-slope propulsion as an individual skill often encountered in the user's environment. According to the WSP, MWUs should keep their weight off the caster wheels by leaning backward. To avoid the downward turning tendency, they should push harder and/or more frequently on the downhill handrim. The WSP states that the uphill hand can be solely used for braking in some cases.

Considering these recommendations and the lack of studies evaluating the propulsion across a cross-slope, we aimed to identify the adaptations of the MWUs to the presence of a crossslope and then to discuss them in light of the WSP recommendations. A secondary objective was to evaluate the potential additional risk of shoulder injury associated with cross-slope propulsion.

METHODS

Study design

This is a cross-sectional exploratory study. The research protocol was approved by the Sainte-Justine University Hospital Center Ethics Committee (MP-21-2020-2533). All participants provided their written informed consent before the experimentation.

Participants

A sample of convenience was recruited from the Institut universitaire sur la réadaptation en déficience physique de Montréal. Interested individuals were included in this study if they had been using a manual wheelchair for mobility for 6 months to 8 years and were 20-to 40-yearold active MWUs with no neurodegenerative disorders. We confirmed that the participants had no history of shoulder pathology based on their self-report.

Experimental Setup

To track the trunk and the right upper limb, 38 reflective skin markers were used. Sixteen additional markers were placed on the wheelchair frame (6) and rear wheels (2x5). Kinematics were recorded using a motion capture system (Vicon Motion Systems, Oxford, UK) at 250 Hz. Muscle activity was recorded using bipolar surface EMG electrodes (Delsys, Natick, USA) for 11 muscles: deltoids (anterior, median, and posterior), trapezius (upper, middle, and lower), serratus anterior, latissimus dorsi, pectoralis major (clavicular head), biceps, and triceps (long and lateral). Additionally, for consenting participants, indwelling electrodes were inserted into the supraspinatus, infraspinatus, and subscapularis. EMG signals were collected at 2000 Hz. Each participant's wheelchair was fitted with a Smart wheel (Three Rivers Holdings LLC, Mesa, USA) on its right side that measured handrim kinetics at 240 Hz.

Procedure

The participants first propelled their wheelchair at a comfortable speed over a level surface. A wooden platform (4.5 m length, 1.2 m width) was used to simulate a 5° cross-slope. The cross-slope angle was chosen based on the set-up used in the Wheelchair Skills Test. [START_REF] Kirby | Wheelchair Skills Test Version 5[END_REF] After getting familiar with the set-up, the participants propelled themselves at a comfortable speed in both directions (right arm on the downhill and uphill of the slope, Fig. 1). They were instructed to maintain a straight-line trajectory. Each condition was repeated four times.

-Insert Figure 1 here -

Additional measures

The details of participants' wheelchairs were logged, as well as the horizontal position of the acromion with respect to the rear wheel axle. Shoulder pain was evaluated using the WUSPI [START_REF] Curtis | Reliability and validity of the Wheelchair User's Shoulder Pain Index (WUSPI)[END_REF] . Neer's [START_REF] Neer | The shoulder in sports[END_REF] and Hawkins-Kennedy 18 tests were used to check for impingement.

Rotator cuff tears were assessed using an external rotation lag sign [START_REF] Hertel | Lag signs in the diagnosis of rotator cuff rupture[END_REF] and a Napoleon test [START_REF] Gerber | Isolated rupture of the subscapularis tendon[END_REF] .

Additionally, participants completed a Wheelchair Skills Test (WST) [START_REF] Kirby | Wheelchair Skills Test Version 5[END_REF] . By interview and selfreported forms, we recorded the daily activities of participants, their major mobility limitations when using a wheelchair, as well as their opinion about cross-slope propulsion.

Analysis

The EMG raw signals were first processed using a fourth-order bandpass filter (20 to 500 Hz and 600 Hz for intra-muscular), then full wave rectified. The envelopes were extracted using a 4 Hz low-pass filter, then normalized to recorded maximal voluntary contractions (MVC).

During the MVCs, the participants increased their activation for 3 seconds before maintaining their MVC for 3 seconds under verbal encouragement. The MVCs were 1 min apart unless the participant required more rest time. For all participants, trunk stability was ensured throughout the MVC. The MVC trials 21 are presented in Table 1. Muscles have been grouped into three categories. The anterior group muscle included the anterior deltoid, the clavicular head of the pectoralis major as well as the serratus anterior. The posterior group consisted of the posterior deltoid, the latissimus dorsi and the trapezius (upper, middle and lower heads).

Finally, the rotator cuff encompassed the infraspinatus, supraspinatus and subscapularis.

-Insert Table 1 here -An Opensim 22 generic model [START_REF] Martinez | Sex differences in upper limb musculoskeletal biomechanics during a lifting task[END_REF] of the trunk and shoulder was scaled to each participant's anthropometry, then supplemented with a basic wheelchair model (frame and two rear wheels). The glenohumeral joint was defined as a ball and socket joint. The joints and wheelchair kinematics were calculated using a least-square inverse kinematics approach [START_REF] Delp | OpenSim: Open-Source Software to Create and Analyze Dynamic Simulations of Movement[END_REF] .

The wheelchair trajectory over the cross-slope was calculated to evaluate how well the participants complied with the instructions. The wheelchair velocity was calculated. The handrim kinetics were first corrected to cancel inaccuracies related to the wheel camber [START_REF] Chénier | Wheelchair pushrim kinetics measurement: A method to cancel inaccuracies due to pushrim weight and wheel camber[END_REF] . The tangential contact forces at the handrim were reported for the downhill and uphill cross-slope propulsion. The trunk orientation with respect to the wheelchair was calculated, as well as the distance between the wheelchair axle and the trunk and upper right-limb center of mass (𝑑 𝐶𝑂𝑀-𝑎𝑥𝑙𝑒 ). This distance represents the wheelchair rotation moment arm around the instantaneous axis of rotation orthogonal to the cross-slope plane. Thus, the bigger the distance, the larger the wheelchair downslope rotation tendency. To account for the difference in wheelchair settings between the participants, the 𝑑 𝐶𝑂𝑀-𝑎𝑥𝑙𝑒 at the resting position was subtracted from the 𝑑 𝐶𝑂𝑀-𝑎𝑥𝑙𝑒 during propulsion. The glenohumeral flexion, abduction and axial rotation were also reported, as well as the elbow flexion. All were compared to their mean values over a level surface. All data analysis was done using Pandas [START_REF]The pandas development team[END_REF] and Numpy [START_REF] Harris | Array programming with NumPy[END_REF] Python libraries.

RESULTS

Nine MWUs participated in this study (Table 2). All participants had spinal cord injuries of grade A or B in the ASIA scale. [START_REF] Kirshblum | International standards for neurological classification of spinal cord injury (Revised 2011)[END_REF] All participants tested negative for impingement and rotator cuff tears. They reported low levels of shoulder pain (<40%) present mostly during wheelchair mobility. All tetraplegic participants used push cuffs. When asked in general about obstacles in urban mobility, none of the participants explicitly spoke about cross-slopes.

However, when asked specifically about cross-slopes, P1, P4 and P6 mentioned no particular difficulty, while the others found the skill demanding, with P3 describing it as uncomfortable.

To mitigate this difficulty, P2 and P7 mentioned switching the sidewalk side, while P9 stated using bike paths when available. Despite the differences in self-reported exertion related to cross-slopes, all participants were able to propel over a cross-slope during the WST, with various success levels: P1, P6, and P9 had trouble maintaining their propulsion direction, P5 had disparate abilities depending on the direction, and P8 applied too much force as his caster wheels lost contact with the ground intermittently. This variability is skill level was also observed for the overall WST score. Finally, while having disparate wheelchair skills training, no participant recalled having training specific for the cross-slope.

-Insert Table 2 here -P5 could not complete all cross-slopes trials as he could not counter the turning tendency sufficiently to keep his wheelchair from rolling off of the cross-slope. Additionally, his latissimus dorsi and serratus anterior EMG data were removed as the sensors lost contact intermittently due to severe lower limb spasms.

Participants propelled their wheelchair slower over a cross-slope (Fig. 2). Their downhill arm had a repeatable push and recovery cycle (Fig. 3). As for their uphill limb, participants' techniques established a spectrum ranging from solely braking the wheel's forward rotation (negative forces as represented by P9) to alternating pushing and breaking the wheel (positive and negative forces as shown for P8).

-Insert Figure 2 here --Insert Figure 3 here -The increased flexion of the trunk increased the 𝑑 𝐶𝑂𝑀-𝑎𝑥𝑙𝑒 (Fig. 4-A). Only P1 had a much more extended trunk over the cross-slope. Whereas all participants showed a lateral trunk flexion to keep their trunk vertical, despite some of them not achieving the same amplitude for both sides (P3 and P5, Fig. 4-B). Participants adapted their glenohumeral kinematics differently to the cross-slope, with an overall tendency to increase the glenohumeral extension and adduction over the cross-slope, as well as the internal rotation for the downhill side (Fig. 5). Overall, the EMG of the anterior group increased for the downhill side, while that of the posterior group increased for the uphill one. Nevertheless, the normalized EMG level remained overall below 60% of the maximal EMG (Fig. 6).

-Insert Figure 4 here --Insert Figure 5 here --Insert Figure 6 here -

DISCUSSION

This exploratory study aimed to identify the adaptations of MWUs to a 5° cross-slope propulsion. Indeed, while few studies investigated cross-slope propulsion, none have investigated explicitly the trunk and uphill arm adaptations. These adaptations are discussed considering the WSP recommendations in the following sections.

Handrim kinetics and wheelchair trajectory

Similar to Sauret et al., [START_REF] Sauret | Handrim mechanical power during wheelchair propulsion on level and cross-slope surfaces: a preliminary study[END_REF] the downhill wheel led the propulsion for all participants and the uphill wheel was modulated to regulate the wheelchair trajectory. Due to the downward turning tendency of the wheelchair, a user must push harder or more often on the downhill wheel. [START_REF] Kirby | Wheelchair Skills Program Manual Version 5.2[END_REF] Participants' approaches to achieve this asymmetry varied in efficiency. Here efficiency refers to the extent to which the work produced by the uphill arm contributes to the execution of the task (i.e., straight propulsion over the cross-slope). P8's technique was the least efficient, as he kept pushing the uphill wheel at each cycle, only to brake its rotation at the "recovery" phase. P8's higher handrim force input did not translate into a faster propulsion, as his propulsion velocity was comparable to that of the other participants. It is interesting to note that P8 already reported shoulder pain during mobility, after only 18 months of wheelchair use. Additionally, while he scored 80% on the WST because he executed advanced skills, he rarely achieved perfect scores, likely relying on strength rather than technique. This strategy might be related to the lack of formal wheelchair skills training.

P9 used his uphill arm to mainly brake the wheel. This technique is included in the WSP tips. [START_REF] Kirby | Wheelchair Skills Program Manual Version 5.2[END_REF] While it should enable the user to reduce energy expenditure to mitigate in part the difficulty of a cross-slope, [START_REF] Brubaker | Effects of side slope on wheelchair performance[END_REF] it might be problematic for tetraplegics. Indeed, tetraplegics rely on applying a large lateromedial force to compensate for their weaker grip. [START_REF] Dallmeijer | Effectiveness of force application in manual wheelchair propulsion in persons with spinal cord injuries[END_REF] Using the downhill arm solely for propulsion, and the upper one exclusively for braking, while applying high enough lateromedial forces to maintain friction between the push-cuff and the handrim, will increase the load on both limbs. As shoulder stability is compromised at high spinal cord injuries, this load would increase the strain on stabilizing muscles. Tetraplegics' difficulty to generate sufficient lateromedial force was observed when P5's arm slipped off the handrim (supplementary material). This suggests the need to assist tetraplegic users through alternative strategies such as using bike paths when possible, as suggested by participants in this study.

Equipment modification (e.g., assistive wheels or projection handrims) should also be investigated as a potential approach. Indeed, projection handrims would ensure control of the wheel's rotation without the necessity of relying on the lateromedial force.

Other participants had some uphill wheel pushes in the first half of the cross-slope. These symmetric pushes might reduce the energy needed to ensure trunk stability against the load transferred through the downhill arm. [START_REF] Astier | Case study: biomechanical analysis of trunk stability in two modes of propulsion of manual wheelchair during start and stabilized speed[END_REF] Additional research is needed to confirm this assumption, and to evaluate the effect of downhill-push and uphill-brake technique on trunk muscles' activation levels.

Trunk kinematics

Reducing trunk flexion over a cross-slope brought the center of mass closer to the rear wheels' axle, which would reduce the turning tendency of the wheelchair. [START_REF] Brubaker | Effects of side slope on wheelchair performance[END_REF] This is coherent with the WSP recommendations. 14 P1's increased trunk extension might have been particularly beneficial as her axle is set posteriorly to her acromion. As most participants' axle was anterior to their acromion, the turning moment arm was already somewhat smaller, which probably eased most participants straight propulsion despite having a trunk flexion similar or superior to that over a level surface. While trunk movement is suspected to contribute to forward propulsion on a level surface, [START_REF] Van Dijk | ISBS Proceedings Archive[END_REF] the reduced propulsion velocity and the turning tendency over the cross-slope probably limit the benefit of such flexion. Additionally, as MWUs who rely on their trunk flexion to push the handrim increase their flexion with fatigue, [START_REF] Rodgers | Influence of trunk flexion on biomechanics of wheelchair propulsion[END_REF] they would struggle to keep their propulsion direction.

For stability, MWUs with motor trunk impairments rely on their latissimus dorsi to provide postural control. [START_REF] Seelen | Impaired balance control in paraplegic subjects[END_REF][START_REF] Louis | Surface electromyography activity of upper limb muscle during wheelchair propulsion: Influence of wheelchair configuration[END_REF][START_REF] Błaszczyszyn | Wavelet coherence as a measure of trunk stabilizer muscle activation in wheelchair fencers[END_REF] They also sit in a kyphotic pose, [START_REF] Minkel | Seating and Mobility Considerations for People With Spinal Cord Injury[END_REF] that leads to a more protracted scapula. [START_REF] Moroder | The influence of posture and scapulothoracic orientation on the choice of humeral component retrotorsion in reverse total shoulder arthroplasty[END_REF] When rolling over a cross-slope, the participants leaned laterally which would induce a transient "scoliotic spine" curvature. This spine configuration decreases scapular range of motion, [START_REF] Rapp Van Roden | The Impact of Trunk Deformity on Shoulder Mechanics and Implications for Upper Extremity Function after Corrective Spinal Surgery in Idiopathic Scoliosis[END_REF] and increases latissimus dorsi activation on the convex side (downhill slope). [START_REF] Kwok | Evaluation of Myoelectric Activity of Paraspinal Muscles in Adolescents with Idiopathic Scoliosis during Habitual Standing and Sitting[END_REF] As the scapula is already protracted, and its range of motion is limited by the induced spine curvature, the upper limb kinematics during propulsion over a cross-slope would rely mostly on the glenohumeral joint, disturbing the scapulohumeral rhythm and increasing the load on scapulohumeral muscles to maintain glenohumeral stability. Additionally, we expect the latissimus dorsi to particularly develop fatigue much faster over a cross-slope, as it would be involved in the lateral spine adjustment, the postural stability, the humeral kinematics, as well as the expiratory function of tetraplegics [START_REF] Fujiwara | EXPIRATORY FUNCTION IN COMPLETE TETRAPLEGICS: Study of Spirometry, Maximal Expiratory Pressure, and Muscle Activity of Pectoralis Major and Latissimus Dorsi Muscles: 1[END_REF] which also increases over a cross-slope. [START_REF] Brubaker | Effects of side slope on wheelchair performance[END_REF] Wheelchair propulsion generates superiorly directed joint reaction forces. [START_REF] Chang | Relationship Between Shoulder Pain and Joint Reaction Forces and Muscle Moments During 2 Speeds of Wheelchair Propulsion[END_REF] After the latissimus dorsi and teres major, the infraspinatus and subscapularis are some of the major contributors to superior stability [START_REF] Halder | Dynamic contributions to superior shoulder stability[END_REF] , putting them at higher risk of fatigue when over a crossslope. On the other hand, compressive forces on the supraspinatus tendon through the humeral head have already been linked to tendinous injury [START_REF] Cook | Is compressive load a factor in the development of tendinopathy?[END_REF] . Thus, the fatigue of the infraspinatus and subscapularis would increase the risk of the supraspinatus tendon degradation.

As the right latissimus dorsi had a higher activation, for most participants, when on the uphill side, the force necessary for the wheel braking was probably higher than for other functions.

From a glenohumeral stability point of view, this activation is appreciated as it would help depress the humerus, reducing the risk of impingement. Indeed, MWUs should be encouraged to actively retract and depress their scapula when braking the uphill wheel, to reduce the risk of large glenohumeral joint reaction forces induced by the posterior deltoid activation.

Glenohumeral joint and muscle excitations adaptations

While each participant modified their glenohumeral kinematics differently over the crossslope, glenohumeral extension and adduction generally increased and internal rotation increased solely for the downhill side. Little can be found in the literature about the risk of shoulder pathology at low elevation angles. For propulsion over a level surface, the humeral head is anteriorly positioned throughout the propulsion cycle, with a low risk of impingement. [START_REF] Mozingo | Comparison of glenohumeral joint kinematics between manual wheelchair tasks and implications on the subacromial space: A biplane fluoroscopy study[END_REF][START_REF] Mozingo | Shoulder mechanical impingement risk associated with manual wheelchair tasks in individuals with spinal cord injury[END_REF] With the glenohumeral adaptations to the cross-slope, the anterior humeral position might increase. [START_REF] Dal Maso | Glenohumeral translations during range-of-motion movements, activities of daily living, and sports activities in healthy participants[END_REF] The compression applied from the humeral head to the biceps long head and supraspinatus tendon is also likely to increase. [START_REF] Cook | Is compressive load a factor in the development of tendinopathy?[END_REF] This is coherent with the conclusions of Sakai et al. [START_REF] Sakai | Characteristic MRI findings of shoulder, elbow, and wrist joints in wheelchair user[END_REF] from MRI observations, relating supraspinatus tears to an upward axial pressure applied by the humeral head. Additionally, the supraspinatus tendon architecture allows for fascicles, separated by a lubricant, to move independently. [START_REF] Fallon | Functional morphology of the supraspinatus tendon[END_REF] This property is probably important over the cross-slope due to the increased tendon excursion. [START_REF] Fallon | Functional morphology of the supraspinatus tendon[END_REF] Yet, the larger compression would oppose the lubrication effect, increasing friction and, consequently, the shear constraint between the fascicles. Such constraint is coherent with the supraspinatus tears occurring at the intrasubstance or articular regions. [START_REF] Morrow | Detailed Shoulder MRI Findings in Manual Wheelchair Users with Shoulder Pain[END_REF] Thus, cross-slope propulsion might increase the anterior glenohumeral head instability and the load on the supraspinatus tendon.

Similar to previous studies that reported no significant differences between level propulsion and a 2.5% cross-slope for the anterior deltoid, the pectoralis major and the infraspinatus peak and mean muscle activations [START_REF] Holloway | Linking wheelchair kinetics to glenohumeral joint demand during everyday accessibility activities[END_REF][START_REF] Symonds | Identifying key experience-related differences in over-ground manual wheelchair propulsion biomechanics[END_REF] , most participants of this study had comparable EMG levels for the leveled and cross-slope propulsions. This lack of difference can be related to the reduced propulsion velocity. Indeed, as participants were instructed to propel themselves at a comfortable velocity, adjusting their velocity might have limited the need to increase their muscle activation. On the other hand, comparing the muscle activation based on the right limb position over the cross-slope, the anterior and posterior muscle groups tend to increase their activation on the downhill and uphill side, respectively. This evolution is particularly visible for the tetraplegic participants. This might highlight the additional burden associated with the limited control of the hand and elbow extensor muscles, as participants try to initiate or brake the rotation of the wheel. Further studies are needed to understand the interaction between the handrim contact forces and the shoulder activation patterns.

WSP and recommendations

Participants were heterogeneous in terms of wheelchair skills with WST scores going from 41% up to 89%. This heterogeneity might be related to how they acquired their skills, as they had various amounts of training. Those with no to little formal training have probably selftrained to acquire more complex skills. Particularly, for cross-slope propulsion, self-training might explain why the participants adapted differently to the cross-slope. As some techniques seemed to be less susceptible to lead to injury and more efficient than others, successfully crossing the slope seems a far cry from simply passing the skill, highlighting the importance of a training program. Given the differences between paraplegic and tetraplegic participants, the training tips for these two groups might differ. Specifically, for those with paraplegia, the existing training tips found in the WSP should be used to improve biomechanical efficiency of this skill, even if the MWU has the strength to do it in a less efficient manner. Indeed, while all paraplegic participants could successfully do a wheelie, none of them opted to propel their wheelchair in a wheelie position over the cross-slope during the WST, despite this technique being the only WSP tip that eliminates the wheelchair turning tendency. For those with tetraplegia, given the additional strain necessary to push or brake the wheel, we recommend providing education to raise awareness of the importance of trunk and glenohumeral stability as well as the benefits of cross-slopes avoidance when possible to minimize their negative impacts. Moreover, our work highlights the need to explore the impact of positioning devices (e.g., lateral supports) and wheelchair modifications (e.g., power assist wheels, handrim projections) on this skill. Such results would help clinicians to weigh the impact of positioning devices and wheelchair modifications on other wheelchair skills or daily-life activities.

As anterior stability is of concern, exercises to help the MWUs identify glenohumeral and thoracohumeral motions are essential to improve the MWUs proprioception so that they start relying on their scapulothoracic muscles. [START_REF] Cools | Rehabilitation of scapular dyskinesis: from the office worker to the elite overhead athlete[END_REF] Particularly, adducting and depressing the scapula should help reduce the strain on the glenohumeral joint.

When using a wheelie is not an option, care should be taken in choosing the technique to achieve asymmetry between the downhill and uphill sides. While alternating pushing and braking on the uphill handrim is probably the worst approach, clinicians should encourage MWUs to try different ratios of push to brake, or different ratios of uphill to downhill push frequencies to find the optimal configuration that would least disturb their trunk stability, account for their strength limitations, and enable trajectory control.

Study Limitations

This study has some limitations. First, as the curvature of the trunk had an impact on the upper limb, and since the adaptations of the trunk probably did not occur at the same spine level, using more markers to track the spine would have given us better insight into trunk adaptations. Secondly, EMG data could not be used to compare participants and should be considered with care. Indeed, some of the muscles were not fully innervated. Thus, a higher EMG does not necessarily equate with a higher muscle force. Testing muscle groups on an isometric ergometer should complete the adaptation's understanding, although it would make for an even longer and more demanding data collection. Our sample of convenience included nine spinal cord injury patients. Accordingly, our conclusions should be further confirmed with a bigger sample size. Nevertheless, this study's results would help guide hypothesisbased future studies to shorten the duration of the data collection, while increasing the number of trials and participants. This study should be replicated for other populations such as patients with stroke as they typically propel with one arm and one leg. It would also be interesting to investigate the effect of the cross slope on the other components of the handrim contact forces. As such study might help further understand the burden associated with crossslope propulsion on the tetraplegic population. Finally, while we measured handrim forces and kinematics for both uphill and downhill sides, our setup was limited to the right limb. It would be interesting to get a synchronous uphill and downhill analysis using experimental data from both the left and right limbs.

CONCLUSIONS

Through this exploratory study, we presented a broader analysis of wheelchair propulsion over a cross-slope, a functional and common task encountered in daily propulsion. The analysis of this skill was done in light of the current clinical recommendations. Based on the analysis of the participants' uphill arm technique and trunk kinematics, we suggested potential injury mechanisms linked to cross-slope propulsion that warrant further investigation and validation. Our study hinted to the importance of wheelchair skills training programs that incorporate condition-specific guidance on posture, kinematics, and muscle control as crucial elements in cultivating a safe and effective technique. Finally, our results should help raise clinicians' and MWUs' awareness to the importance of exploring wheelchair modifications and positioning devices to reduce potential for injury.

FIGURES AND TABLES

POSITION RESISTANCE

Seated with the arm flexed at 90°, elbow extended, and palm facing down.

Resistance to arm flexion at the elbow.

Seated with the arm abducted at 90°, elbow extended, and palm facing down.

Resistance to arm abduction at the elbow.

Prone with the arm horizontally abducted at 90° and the elbow flexed at 90°. Resistance to horizontal abduction at the elbow Seated with the arm flexed at 125° and the elbow extended.

Resistance to arm flexion at both the elbow and the inferior angle of the scapula.

Seated with the arm abducted at 30° and the elbow flexed at 90°. Resistance to arm adduction at the elbow.

Seated with the arm abducted at 90°, the elbow extended, the palm facing down, and the neck bent to the arm side.

Resistance to arm abduction at both the head and the elbow. Prone with the arm at 90° horizontal abduction, the elbow extended and the thumb pointing to the ceiling.

Resistance to horizontal abduction at the elbow.

Prone with the arm horizontally flexed at 90° in line with the lower trapezius line of action, the elbow extended and the thumb pointing to the ceiling. Resistance to horizontal flexion at the elbow.

Seated with the arm abducted at 90° and the elbow flexed at 90°. Resistance to arm adduction at the elbow.

Seated with the arms flexed at 90° and the elbows flexed so that the palms of both hands are in contact.

Pressing hands together.

Lying on the side with the arm at the side, the elbow extended and the palm facing down.

Resistance to arm abduction at the wrist.

Lying on the side with the arm at the side and the elbow at 90° flexion.

Resistance to external rotation at the wrist.

Prone, with the back of the hand in contact with the upper lumbar spine. Resistance to internal rotation at the wrist.

Seated, the arm at a neutral position and the elbow at 90° flexion with the forearm in supination.

Resistance to elbow flexion at the wrist.

Seated, the arm 90° abducted, the elbow at 90° flexion with the forearm in a neutral. Resistance to elbow extension at the wrist.

Fig. 1 .Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .

 123456 Fig. 1. Conditions of the wheelchair propulsion

  Notes: WST-CS is the WST score for the cross-slope skill. Details about the WST and the WUSPI scoring can be found in supplementary material. MWC training reports the answer of the participants when asked about their training during the interview. No training meant no formal training. Abbreviation: f(emale), m(ale), MWC(manual wheelchair) * The injury level assessment evolved from a C4-5 to C6-7 since the time of injury.

Table 1 :

 1 Participant's position and applied resistance during each MVC. 1

Table 2 .

 2 Demographics and characteristics of the participants P1 to P9 2

	Participant (sex)	Age (years)	Mass (kg)	Months using the MWC	MWC mass (kg)	Injury level	WUSPI [0-150]	WST (%)	WST-CS [0-3]	MWC training
	P1 (f)	29	73	48	12.9	T10	19	60.22	2	No training
	P2 (m)	34	67	72	13.5	C7	19.28	87.10	3	~ 10 h of training
	P3 (m)	28	70	48	22.2	C4-5	40	41.38	2	3-4 h of training
	P4 (m)	40	72	12	14.4	T8	15	68.82	2	Two sessions a week for about a year
	P5 (m)	24	68	41	16.5	C4-5 (C6-7) *	46.25	50.00	1	No training
	P6 (m)	27	65	24	19.1	T3	24.65	75.27	2	Training every day for 3 months
	P7 (m)	38	58	18	16.8	T11-12	62.5	88.89	3	Some average training
	P8 (m)	29	80	18	15.1	T4	29	80.00	2	No training
	P9 (m)	39	89	19	16.9	C6-7	26.25	47.78	2	Few average 15-min training sessions
	Mean ± std	32 ± 5.5	71.3 ± 9.1	33.3 ± 18.9	16.4 ± 2.7		33.8 ± 14.3 66.61 ± 16.61

This research was undertaken thanks, in part, to funding from the Canada First Research Excellence Fund through the TransMedTech Institute. This study was carried out within the Associated International Laboratory EVASYM. The authors acknowledge Frédéric Messier and the Centre intégré universitaire de santé et de services sociaux du Centre-Sud-de-l'Ile-de-Montréal for their help in participant recruitment. no conflict of interest to declare.

Trunk and glenohumeral joint adaptations to manual wheelchair propulsion over a cross-slope: an exploratory study

Appendix A: WST and WUSPI scores

The WST includes multiple skills that are scored individually. Each skill is scored on capacity from 0 to 3. Briefly, a score of three indicates an advanced pass, a score of 2 indicates that the participant meets all evaluation criteria, but still has room for improvement. A score of 1 is given when the subject completes the majority but not all the evaluation criteria of the given skill. Finally, a score of 0 is given when the participant's execution of the skill is incomplete, unsafe, or if the participant is unwilling to try. A skill can be evaluated as "Not Possible" if the wheelchair does not have all parts that would allow for a safe execution, or as "Testing

Error" if the evaluator is incapable of assessing the skill. Finally, the global WST score is evaluated as a percentage of success in individual skills other than those scored as "Not

Possible" or as "Testing Error". The formula is shown below:

∑ 𝑖𝑛𝑑𝑖𝑣𝑖𝑑𝑢𝑎𝑙 𝑠𝑘𝑖𝑙𝑙𝑠 𝑠𝑐𝑜𝑟𝑒𝑠 (𝑛 𝑠𝑘𝑖𝑙𝑙𝑠 -𝑛 𝑠𝑘𝑖𝑙𝑙𝑠 𝑠𝑐𝑜𝑟𝑒𝑑 "𝑁𝑜𝑡 𝑝𝑜𝑠𝑠𝑖𝑏𝑙𝑒" -𝑛 𝑠𝑘𝑖𝑙𝑙𝑠 𝑠𝑐𝑜𝑟𝑒𝑑 "𝑇𝑒𝑠𝑡𝑖𝑛𝑔 𝐸𝑟𝑟𝑜𝑟" ) * 3 * 100

The WUSPI is a 15-item self-report measure of shoulder pain in wheelchair during functional daily activities. The activities are grouped under four categories: transfers, wheelchair mobility, self-care and general activities. For each activity, the participant evaluated their pain between "no pain" (0) to "worst pain ever" (10). The sum of the individual activities' scores is the raw WUSPI score. To account for activities that were not applicable for certain participants, a performance corrected score was obtained by dividing the raw WUSPI score by the number of applicable activities then multiplying it by 15.