Daniele Faenzi
email: daniele.faenzi@u-bourgogne.fr

Marcos Jardim
email: jardim@unicamp.br

Jean Vallès
email: jean.valles@univ-pau.fr

Jean Vallès The

Saito Criterion

Jean Vall Ès

Enrique Arrondo

THE SAITO CRITERION AND SOME AVATARS

Keywords: 2010 Mathematics Subject Classification. AF404, 14J60, 14M10, 32S65 Logarithmic module, freeness and local freeness. Complete intersection

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Let R = k[x 1 , . . . , x n ] be a graded commutative ring over a field k of characteristic zero. The partial derivatives ∂ ∂x i are denoted ∂ i . The R-module Der R of polynomial derivations is a free rank n module:

Der R ≃ ⊕ 1≤i≤n R∂ i .

To a square-free polynomial f ∈ R, one associates its module of tangent derivations (or module of tangent vector fields):

Der(f ) = {δ ∈ Der R | δ(f ) ∈ (f )},
where (f ) is the principal ideal generated by f . The module Der(f ) is a rank n reflexive R-module (see [START_REF] Saito | Theory of logarithmic differential forms and logarithmic vector fields[END_REF]Cor. 1.7]) that could also be defined by the canonical exact sequence

0 -→ Der(f ) -→ Der(R) ≃ R n -→ R/(f )
where the last map is

Der(R) ≃ R n -→ R/(f ), δ → δ(f ) mod (f ).
Remark 1.1. When V (f ) is a normal crossing divisor (NCD for short) then Der(f ) is locally free (see [START_REF] Deligne | Équations différentielles à points singuliers réguliers[END_REF] for instance). However, the converse is not true.

Remark 1.2. The map Der(R) ≃ R n -→ R/(f ) is the composition of R n -→ R sending δ on δ(f ) and the quotient map R -→ R/(f ) sending a polynomial on its reduction modulo f . The image of the first map is the Jacobian ideal J ⊂ R generated by the partial derivatives ∂ i f . Its reduction modulo (f ) is J+(f ) (f ) ; when f is homogeneous, J + (f ) = J. It defines a subscheme in V (f ) supported by the singular locus of V (f ). This locus lives in codimension > 1 since f is square free.

Remark 1.3. If f is homogeneous then δ E := i x i ∂ i ∈ Der(f ).
The derivation δ E is known as Euler derivation. This induces a splitting

Der(f ) = Rδ E ⊕ Der 0 (f ) where Der 0 (f ) = {δ ∈ Der(R) | δ(f ) = 0}; to be precise, if f is homogeneous and δ ∈ Der(f ) then δ -1 deg(f ) δ(f ) f δ E ∈ Der 0 (f ). The module Der 0 (f ) is the Syzygy module of ∇f = (∂ 1 f, . . . , ∂ n f ), 0 -→ Der 0 (f ) -→ R n ∇f -→ R
Remark 1.4. Der(f ) being reflexive, it must be locally free on the affine plane A 2 (resp. P 2 ).

Remark 1.5. When V (f ) is an hyperplane arrangement in P n , that is V (f ) = ∪ i H i is a union of N + 1 distinct hyperplanes in P n then T f = Der(f ) is the dual of a Steiner sheaf, more precisely it arises as the kernel of a matrix of linear forms that could be described by the linear relations between the H i (see [START_REF] Dolgachev | Arrangements of hyperplanes and vector bundles on P n[END_REF], [START_REF] Faenzi | Logarithmic bundles and line arrangements, an approach via the standard construction[END_REF], [START_REF] Arrondo | Schwarzenberger bundles of arbitrary rank on the projective space[END_REF], [START_REF] Ancona | Unstable hyperplanes for Steiner bundles and multidimensional matrices[END_REF], etc.):

0 -→ T f -→ O(-1) N -1 -→ O N -n-1 .
As first pointed out by Saito in [START_REF] Saito | Theory of logarithmic differential forms and logarithmic vector fields[END_REF], there are divisors V (f ) not NCD such that Der(f ) is a free R-module (meaning that Der(f ) ≃ R n ). In the same paper Saito gave a determinantal characterization of freeness, commonly called Saito criterion. The aim of this note is to prove again this criterion (see Theorem 3.3) with basic algebraic arguments in order to generalize it to non reduced divisor (see Theorem 4.1) and overall to finite families of polynomials (see Theorem 6.1). All these avatars of Saito criterion are direct consequences of a elementary statement about reflexive module which is explained in the next section 2.

Except when it is specified, the polynomials considered in this text are not necessarily homogeneous.

Algebraic prelude: key proposition

We begin by a elementary proposition concerning reflexive modules. This proposition will be the key argument in the next proofs of Saito criterion and its generalization.

Proposition 2.1. Let M be a rank s ≤ n reflexive R-module and Φ a monomorphism M ֒→ R n . Let N = 1≤i≤s Re i be a free sub-module of M and ι : N ֒→ M be the inclusion map. Then

s (e 1 , . . . , e s ) = s Φ • ι = det(ι) s Φ.
In particular, the following statements are equivalent:

(1)

M = N (2) det(ι) ∈ k * (3) V ( s (Φ • ι)) = V ( s Φ)
Proof. This is a consequence of the commutative square:

0 N ι N (e 1 ,...,es) 0 / / M Φ / / R n (1)
Indeed, N and M being reflexive modules of rank s, s N ≃ s N ∨∨ and s M ≃ s M ∨∨ are reflexive module of rank one that is free modules and there is a commutative square

0 s N det(ι) s N s (e 1 ,...,es) 0 / / s M s Φ / / s R n (2)
showing that

s (e 1 , . . . , e s ) = s (Φ • ι) = det(ι) s Φ.
This proves our proposition. Let us point out that if s = n then det(Φ) | det(e 1 . . . , e s ).

Free divisors: Saito criterion

Definition 3.1. The module Der(f ) is free of rank n if and only if there exist n derivations

δ i = 1≤j≤n P ij ∂ j ∈ Der(f ) such that Der(f ) = ⊕ i Rδ i .
Let us first point out the following crucial fact (see [Sai80, 1.5]):

Lemma 3.2. Let δ i ∈ Der(f ) be n independent derivations tangent to V (f ). Then, f | det(δ 1 , . . . , δ n ).
Proof. Since δ i are n independent derivations tangent to V (f ) the R-module ⊕Rδ i is a free sub-module of rank n of Der(f ). Calling N = ⊕Rδ i and M = Der(f ) we have, according to Remark 1.2, det(Φ) = f and according to Proposition 2.1,

det(Φ) = f | det(Ψ) = det(δ 1 , . . . , δ n )
As direct consequence we get the classical Saito criterion (see [Sai80, Theorem (ii), page 270]): Proof. Let U = ⊕ j=1,...,n Rδ j be a free sub-module of Der(f ). We get a commutative diagram:

Theorem 3.3. Der(f ) = ⊕ i Rδ i if and only if det(δ 1 , . . . , δ n ) = uf, u ∈ k * .
0 0 / / U / / R n / / K / / 0 / / Der(f ) / / R n / / J+(f ) (f ) / / 0 Der(f )/U 0 0 (3)
where the singular locus of the cokernel K is defined by

V ( n U ) = V (det(δ i )).
Since Der(f ) is reflexive, U is free and rk(Der(f )) = rk(U ), then Der(f )/U = 0 when Der(f ) = U or Der(f )/U is a torsion sheaf when Der(f ) = U.

-If U = Der(f ) then it is is free and K = J+(f ) (f ) by the snake lemma. This proves that the singular locus V ( n U ) of K is exactly V (f ). Actually, this proves more. It proves that the Jacobian ideal J modulo f coincide with the first Fitting ideal defined by V ( n-1 U ).

-If Der(f ) = U then Der(f )/U is a torsion sheaf. By the snake lemma we get

0 -→ Der(f )/U -→ K -→ J + (f ) (f ) -→ 0
proving that the singular locus of K defined by V (det(δ i )) differs from V (f ).

Multiderivations

The previous material can be extend to non reduced divisors and their associated modules of multiderivations. Indeed a multiple polynomial

f = f m 1 1 • • • f mr r
where f i are square free polynomials of degree r i ≥ 1 being given one can define the following sub-module of Der(R) (see [START_REF] Günter | Multiarrangements of hyperplanes and their freeness, Singularities[END_REF]):

Der m (f ) := {δ ∈ Der R | δ(f i ) ∈ (f m i i ) 1 ≤ i ≤ r}.
This R module which is also reflexive, can be defined in by the canonical exact sequence

0 -→ Der m (f ) -→ Der(R) ≃ R n -→ R/(f )
where the last map is

Der(R) ≃ R n -→ R/(f ), δ → δ( i f i ) mod ( i f m i i )).
In particular this shows that its rank is n. Moreover, there is a natural inclusion

Der m (f ) ⊂ Der( i f i ).
Saito's criterion is still pertinent for this module, as this is noted in [START_REF] Günter | Multiarrangements of hyperplanes and their freeness, Singularities[END_REF]:

Theorem 4.1. Der m (f ) = ⊕ i Rδ i if and only if det[δ 1 , . . . , δ n ] = uf, u ∈ k * .
Proof. Same proof than Theorem 3.3.

4.1. Points in P 1 . Assume now that the polynomials f i ∈ k[x, y] are linear forms (i.e. homogeneous of degree 1). Then V (f ) consists in r distinct points p 1 , . . . , p r with multiplicity m 1 , . . . , m r .

Lemma 4.2.

Der(f 1 • • • f r ) = R[-1] ⊕ R[1 -r].
Proof. The Euler derivation and the Koszul derivation

∂ y f ∂ x -∂ x f ∂ y define a free sub- module in Der(f 1 • • • f r ). By Theorem 3.3 this sub-module is the whole Der(f 1 • • • f r ). Since Der m (f ) is reflexive it must be free on R = k[x, y]. Then there are two integers 0 ≤ a ≤ b such that Der m (f ) = R[-a] ⊕ R[-b] and a + b = i m i . Since Der m (f ) ⊂ Der(f 1 • • • f r ) we have an exact sequence: 0 -→ R[-a] ⊕ R[-b] -→ R[-1] ⊕ R[1 -r] -→ ⊕ i R/(f m i -1 i )

Examples of applications

We give here two kinds of applications, some in the homogeneous (or projective) case and some in the non homogeneous case (or affine). The first case is in some sense easier since the degree of the derivations are or are not the good ones, separating the free from the non free divisors. In the affine case, the criterion appears to be stronger since the degrees are not enough to distinguish which set of derivations is a basis or not.

Example 5.1. The arrangement of 3n + 3 lines (the reflection arrangement denoted by

A 3 3 (n) in [OT92]) on P 2 defined by the equation f (x, y, z) := xyz(x n -y n )(x n -z n )(y n -z n ) = 0
is free with exponents (n + 1, 2n + 1). Indeed the derivations

δ 1 = x n+1 ∂ x + y n+1 ∂ y + z n+1 ∂ z and δ 2 = x 2n+1 ∂ x + y 2n+1 ∂ y + z 2n+1 ∂ z belong to Der(f ) and det   x x n+1 x 2n+1 y y n+1 y 2n+1 z z n+1 z 2n+1   = xyzdet   1 x n x 2n 1 y n y 2n 1 z n z 2n   = xyz(x n -y n )(x n -z n )(y n -z n ).
Let us give two examples now of smooth conics in the affine plane having a basis consisting in polynomials of degree 2 for the first one and of degree one for the other.

Example 5.2. This exemple is due to Yohan Genzmer. He shows that even if Der(f ) contains a derivation of minimal degree this derivation is not necessarily a generator of Der(f ). Let f (x, y) = x 2 + y 2 + x. Then the derivation 2y∂ x + (2x + 1)∂ y ∈ Der(f ) in particular it belongs to Der 0 (f ) but a basis of Der(f ) is given by

δ 1 = 2y(2x + 1)∂ x + (4y 2 -1)∂ y δ 2 = x(x + 1)∂ x + ( y 2 + yx)∂ y Indeed δ 1 (f ) = 2yf , δ 2 (f ) = (2x + 1)f and det(δ 1 , δ 2 ) = f. Actually the curve V (f ) ⊂ A 2 is free but the curve V (x 2 + y 2 + xz) ⊂ P 2 is not free.
Example 5.3. The curve C ∪ L = {z(x 2 + yz) = 0} union of a smooth conic and one tangent line in P 2 is free with exponents (1, 1). Indeed two derivations in Der 0 (C ∪ L) are z∂ x -2x∂ y and x∂ x + 4y∂ y -2z∂ z . and det

  x z x y -2x 4y z 0 -2z   = 6z(x 2 + yz).
Removing the line L = {z = 0} we get

0 -→ T C∪L = O P 2 (-1) 2 -→ T C -→ O L (α) -→ 0,
for some integer α. Then on the affine space P 2 \ L, the curve becomes {y + x 2 = 0} which is free with two derivations of degree 1:

δ 1 = ∂ x -2x∂ y δ 2 = x 2 ∂ x + y∂ y verifying det(δ 1 , δ 2 ) = y + x 2 .

Saito criterion for several polynomials

Let S = V (f 1 , . . . , f s ) be a locally complete intersection of 1 < s ≤ n -1 reduced hypersurfaces. We denote by I S = (f 1 , . . . , f s ) the ideal defining S. The module of logarithmic derivations associated to S, say Der(S), consists in derivations δ such that δ(g) ∈ I S for all g ∈ Der(S). This is a rank n module; one can verify easily that it contains Der( i f i ).

In this section we are interested in the module:

Der(f 1 , . . . , f s ) = {δ ∈ Der(R) = R n | f j δ(f i ) = f i δ(f j ), 1 ≤ i < j ≤ s}.
We remark first that this module depends on the choice of the polynomials generating I S .

Theorem 6.1. Let f 1 , . . . , f s ∈ R be 1 < s ≤ n -1 algebraically independent polynomials. Then, (1) Der(f 1 , . . . , f s ) is a reflexive module of rank n -s + 1.

(2) Let δ 1 , . . . , δ n+1-s be independent derivations in Der(f 1 , . . . , f s ). Then,

Der(f 1 , . . . , f s ) = ⊕ j=1,...,n+1-s Rδ j ⇔ codim k n V n+1-s [δ 1 , . . . , δ n+1-s ] > 1.
Proof. Let us consider the R-module map:

DerR = R n -→ R s-1 , µ → (f i µ(f 1 ) -f 1 µ(f i )) 2≤i≤s .
The image of this map, let say K, being a sub-module of a free module is torsion free. Since the f i are algebraically independent this map is of maximal rank. Being a map between two free modules, its kernel is reflexive. Moreover, the kernel of this map is the module Der(f 1 , . . . , f s ) since it coincides with the kernel of the map

DerR = R n -→ R ( s 2 ) , µ → (f i µ(f j ) -f j µ(f i ))
1≤i<j≤s . This proves (1).

Considering now the free sub-module

U = ⊕ j=1,••• ,n+1-s Rδ j of Der(f 1 , . . . , f s ) we obtain a commutative diagram: 0 0 / / U Ψ / / ι R n / / K ′ / / 0 0 / / Der(f 1 , . . . , f s ) Φ / / R n / / K / / 0 Der(f 1 , . . . , f s )/U 0 (4) Since K is torsion free this implies that codim k n V n+1-s Φ > 1.
Moreover, by Proposition 2.1, n-s+1 Ψ = det(ι) n-s+1 Φ proving that U = Der(f 1 , . . . , f s ) ⇔ det(ι) ∈ k * ⇔ codim k n V n+1-s Ψ > 1.

Example 6.2. Let us give an example in A 3 . Let f (x, y, z) = x 2 -y 2 and g(x, y, z) = z 3 +1. The image of the map R 3 -→ R, µ → f µ(g) -gµ(f ) is the ideal (x(z 3 + 1), y(z 3 + 1), z 2 (x 2 -y 2 )). The derivations

δ 1 = y∂ x + x∂ y δ 2 = xz 2 ∂ x + yz 2 ∂ y -(z 3 -1)∂ z
belong to Der(f, g) (actually δ i (f ) = δ i (g) = 0). Since codimV ( 2 (δ 1 , δ 2 )) > 1, the module Der(f, g) is free, more precisely, Der(f, g) = Rδ 1 ⊕ Rδ 2 .

Example 6.3. The previous example comes from the simple example in P 3 given by the two equations: f = x 2 0 -x 2 1 , g = x 3 2 + x 3 3 . Since the syzygies of the Jacobian ideal are generated by (x 1 , x 0 , 0, 0) and (0, 0, x 2 3 , -x 2 2 ) the sheaf T (f,g) = O(-1) ⊕ O(-2).

If we take, instead of g, another cubic in the ideal (f, g) for instance g ′ = g + x 0 f then T (f,g ′ ) is no more free because there is no syzygy of degree 1. If we take g" = g + x 2 f then it is still free and verifies also T (f,g) = O(-1) ⊕ O(-2).

3. 1 .

 1 First proof of Theorem 3.3. Proof. Denoting N = ⊕Rδ i and M = Der(f ) we have, according to Proposition 2.1, M = N if and only if det(ι) = u ∈ k * . 3.2. Second proof of Theorem 3.3.

D.F. partially supported by FanoHK ANR-20-CE40-0023, SupToPhAG/EIPHI ANR-17-EURE-0002, Région Bourgogne-Franche-Comté. M.J. is supported by the CNPQ grant number 305601/2022-9 and the FAPESP Thematic Project number 2018/21391-1. All authors partially supported by Bridges ANR-21-CE40-0017 and CAPES/COFECUB project Moduli spaces in algebraic geometry and applications, Capes reference number 88887.191919/2018-00.


Mathias Samson

Anthony Bernard

Tony Monnet

Patrick Lacouture

Laurent David

Unsteady computational fluid dynamics in front crawl swimming

Keywords: URANS method, CFD, propulsive forces, k-ω SST turbulence model, VOF method

  L'archive ouverte pluridisciplinaire

Introduction

In competitive swimming, the hand and forearm play a major role in the propulsion of the swimmers [START_REF] Maglischo | Swimming Fastest[END_REF][START_REF] Takagi | Measurement of propulsion by the hand during competitive swimming[END_REF]. Moreover [START_REF] Berger | Hydrodynamic drag and lift forces on human hand/arm models[END_REF], [START_REF] Toussaint | Pumped-up propulsion" during front crawl swimming[END_REF] and [START_REF] Kudo | The effect of unsteady flow due to acceleration on hydrodynamic forces acting on the hand in swimming[END_REF] highlight the leading role of the hand relative to the forearm (67 % of the hand-forearm propulsive force according Schleihauf et al., 1983). The propulsive forces are created by the flow surrounding the swimmer's arms. The hand motions in front crawl swimming are complex: it is three-dimensional with sudden changes in velocities and angles of attack [START_REF] Samson | Kinematic hand parameters in front crawl at different paces of swimming[END_REF].

Thus, the induced flow is highly non-linear and turbulent containing complex unsteady vortex dynamics [START_REF] Toussaint | Pumped-up propulsion" during front crawl swimming[END_REF][START_REF] Polidori | Skin-friction drag analysis from the forced convection modeling in simplified underwater swimming[END_REF]Zaïdi et al. 2010;Takagi et al. 2015). It is therefore extremely difficult to study, both experimentally and numerically. Nevertheless, progress in recent years, not only in the development of numerical methods and turbulence modelling, but also in computing power, currently allows to simulate more complex flows. CFD is seen as an interesting tool because it allows to simulate, the flow around a given geometry, and the external forces applied to the solid parts. [START_REF] Bixler | Computational Fluid Dynamics: An Analytical Tool for the 21st Century Swimming Scientist[END_REF] have shown that studies conducted in CFD in quasi-static conditions are insufficient to model a complex flow around a swimmer. Indeed, these methods take into account the average of the flow field, without taking into account the unsteady effects, engendered by the fast instantaneous velocity variations [START_REF] Sato | Computational Fluid Dynamics Analysis of Hydrodynamic Force Acting on a Swimmer's Hand in a Swimming Competition[END_REF][START_REF] Rouboa | The effect of swimmer's hand/forearm acceleration on propulsive forces generation using computational fluid dynamics[END_REF]). These effects, however, play a major role in the overall dynamics of the flow. The first unsteady effect highlighted in swimming was the mass effect added [START_REF] Pai | A hydrodynamic study of the oscillation motion in swimming[END_REF]. Ellington in 1984 defines it as is the reaction force due to the acceleration and deceleration of the fluid particles around the moving solid. A second effect (bound vortex), which is the circulation around a body following a vortex shedding, plays an important role in swimming, especially during the phases of change of direction of the hand [START_REF] Takagi | Unsteady hydrodynamic forces acting on a robotic arm and its flow field: Application to the crawl stroke[END_REF]). These two effects favorably contribute to the thrust.

Recently, unsteady studies have been proposed to analyse the flow around the swimmers. [START_REF] Lecrivain | Effect of body roll amplitude and arm rotation speed on propulsion of arm amputee swimmers[END_REF] have numerically investigated the effect of body roll amplitude and of upper arm rotation speed on the propulsion of an arm amputee swimmer. [START_REF] Keys | A Full Body Computational Fluid Dynamic Analysis of the Freestyle Stroke of a Previous Sprint Freestyle World Record Holder[END_REF], from a smoothed particle hydrodynamics (SPH) method, examined propulsive and drag forces experienced across the body during full-body freestyle swimming. Both used the morphing mesh method to model the joint segment, with the k-ε turbulence model. Moreover, [START_REF] Sato | Computational Fluid Dynamics Analysis of Hydrodynamic Force Acting on a Swimmer's Hand in a Swimming Competition[END_REF] simulated, by a moving mesh method, a flow surrounding a hand but without turbulence model or free surface.

Turbulence models and the local Reynolds number are very important in the simulation because the role played by the boundary layer in the creation of turbulence. Turbulence models have been improved a lot in recent years, in particular by taking more account of near wall phenomena (Manceau 2010). So, the turbulence models used in the research of swimming have progressed from a k-ε turbulence model with Bixler (2002), to a k-ω model more recently with Popa et al. (2014) and [START_REF] Arfaoui | Swimming drafting simulation using a k-omega turbulence model[END_REF]. These two eddy-viscosity models are the most popular in CFD studies in this domain. [START_REF] Banks | An analysis of a swimmer's passive wave resistance using experimental data and CFD simulations[END_REF] used the Shear Stress Transport (SST) k-ω turbulence model [START_REF] Menter | Two-Equation Eddy-Viscosity Turbulence Models for Engineering applications[END_REF] to measure the passive drag around a swimmer. This model better reflects the physics than k-ε model, because it combines the qualities of the high Reynolds k-ε model (away from walls), and the low Reynolds k-ω model (near the walls).

In this context, an unsteady CFD study is presented that applies the URANS method (Unsteady Reynolds Averaged Navier-Stokes). This method is widely used in aeronautics, and although it remains a subject of active research, it has been validated in many canonical unsteady flows: wake, downwardstep, and compressible co-axial jets [START_REF] Fadai-Ghotbi | Modélisation de la turbulence en situation instationnaire par approches URANS et hybride RANS-LES. Prise en compte des effets de paroi par pondération elliptique[END_REF]. To better approximate the physics of the flow, we propose a free surface modeled by a Volume of Fluid method (VOF), and we use the k-ω SST turbulence model. We employ a polyhedral finite-volume Navier-Stokes flow solver, coupled with of an overset mesh method. The main purpose was dual: first to analyse the evolution of the hand-forearm propulsive forces, and explain this relative to the arm kinematics parameters. Second, to validate this model by comparaison with others results.

Methods

Numerical method

URANS method

The flow around a swimmer is turbulent and unsteady with high Reynolds number, thus a full Direct Numerical Simulation (DNS) is not possible with today's computing power. A method to resolve this shortcoming is to use the Reynolds averaged Navier Stokes equations, closed by a turbulence model, under the unsteady form (URANS equations). These equations are obtained by introducing the Reynolds decomposition, which takes into consideration that, in turbulent flows, each instantaneous variable is the sum of a mean component and a fluctuating component (equations ( 1), (2)):

) , ( + > ) , ( =< ) , ( ' t x u t x U t x U i i i (1) ∑ 0 = ∞ → ) + , ( 1 + 1 lim >= ) , ( < N n i N i nT t x U N t x U (2) ) t , x ( i U is the instantaneous velocity of the hand component in the i-direction (m/s),   ) t , x ( i U is the phase average of velocity component in the i-direction (m/s), ) t , x ( ' i u
is the fluctuating velocity component in the i direction (m/s). Equation ( 1) is obtained from the double decomposition [START_REF] Cantwell | An experimental study of entrainment and transport in the turbulent near wake of a circular cylinder[END_REF]. In the equation 2, T is the period of the cyclic phenomenon, i.e. from the entry of the hand into the water to the exit.

Then, the time averaging of the instantaneous equations leads to the following system of averaged equations ((3) and ( 4)). These equations are formally identical to the RANS equations for a statistically unsteady flow, taking into account the unsteady term t ∂

U ∂ i   . 0 i x i U =     (3) (4)
This URANS decomposition can be interpreted as the application of an implicit spatio-temporal filter at instantaneous variables, this filter being employed by the introduction of a turbulence model in the

Turbulence model

The k-ω SST turbulence model was used for closing the URANS equations. This is a two-equation eddyviscosity model, in which k is usually interpreted as the turbulence kinetic energy and ω is the rate of dissipation of the eddies. k-ω SST is a "hybrid" of two models: the k-ω model [START_REF] Wilcox | Reassessment of the scale-determining equation for advanced turbulence models[END_REF] which is efficient in the near-wall region (but not applicable to external aerodynamics) and the standard k-ε model [START_REF] Launder | The numerical computation of turbulent flows[END_REF] which gives good results far from the wall. The transition between the two models is made by quite complex functions which depend on the distance to the wall. This model has been preferred to other low-Reynolds models (as the low-Reynolds k-ε model for instance) because k-ω SST is integrable down to the wall, without introducing a damping function, which allows to take more account of the geometry of the solid. Moreover the equation of ω is independent of that of k, which is a factor of numerical robustness.

These parameters are defined from two transport equations:

(5) (6) (kinematic eddy viscosity) is calculated as a function of k and ω by the following equation:

The use of a k-ω SST formulation in the inner parts of the boundary layer makes the model directly usable all the way down to the wall through the viscous sub-layer. It requires a higher resolution mesh to solid boundaries (see below) with the first calculation point in the viscous sublayer.

Flow solver

The numerical simulations have been carried out with the STAR-CCM+ ® software. The finite volume method is applied on a hybrid mesh made of structured mesh (near-wall) and unstructured mesh. The computation is made from a computer with 736 processors at 2.7 GHz and 1.47 TB of RAM. The linkage between the momentum and continuity equations is achieved with a predictor-corrector approach, with a standard SIMPLE algorithm. The advancement in time is done by an implicit first-order Euler scheme.

This physical time-step was 2.5.10 -3 s (with 30 inner iterations) and was chosen to obtain a Courant-Friedrichs-Lewy number (𝑈∆ 𝑥 ∆ ⁄ 𝑡) close to 1 when U is maximal.

Free surface

The free surface is simulated with multi-phase liquid and air of constant density from a Volume of Fluid method (VOF). According to [START_REF] Dabnichki | Unsteady fluid mechanics effects in water based human locomotion[END_REF], the presence of the free surface should affect the force generation especially in the initial and final phases when the arm is close to the surface.

Overset grid method

This method consists to move a mesh (overset) into a fixed mesh (background). The overset grid method enables to efficiently simulate a movement of solid body (here the hand-forearm) in a fluid. The volumes control domain consisted of two regions (Figure 1): the background (parallelepiped 5m x 3m x 1m, containing 280,000 cells), fixed, and an overset (sphere of 1m in diameter containing 820,000 cells), attached at the hand-forearm and moving with it with six degrees of freedom (three in rotation and three in translation). The mesh consisted of polyhedral and pyramidal cells. It was refined and structured at the near wall of the hand-forearm, using 10 prism layers, stretching at 1.1, which the size of the first cell is equal to 10 -4 m to obtain y + ≈ 1. y + being a dimensionless number equal to y𝑢 𝜏 𝜈 ⁄ , with 𝑢 𝜏 the friction velocity ( 𝑢 𝜏 = √𝜏 𝑝 𝜌 ⁄ , 𝜏 𝑝 : wall shear stress, 𝜌: fluid density).

Figure 1 The mesh is also refined at the free surface. Beyond the near wall areas, the mesh is unstructured in the overset volume. The background mesh is structured, with average mesh size equal to 0.07 m. The independence of the calculated physics based on the numerical conditions was tested with respect to the computational time-step, the number of iterations, and the domain and mesh size. The interface between the two regions was then calculated by interpolation between active cells near each region.

Boundary conditions

Two inlet boundaries were defined (Figure 1). A first, vertical, composed of an air / water in variable composition (managed by a field function), the input velocity is equal to 0 m.s -1 for reproducing the condition of a fluid at rest. The second (top surface) consists only of air. An outlet boundary (in pressure condition) is composed of an air / water varying in composition to ensure the conservation equations.

Both sides of the volume are defined as a symmetry plane. The surface of the arm and the bottom of the volume are defined as a wall to obtain no-slip conditions. The free surface is set at a height of 2 m above the pool bottom.

Kinematic parameters of the simulation

Geometry of the hand-forearm surface

The simulation consisted of reproducing the trajectory of a hand-forearm of a swimmer in sprint pace, from experimental kinematics measurements. The hand and forearm segment (with the wrist locked) of an expert swimmer (level: 90.2 % of world record, height: 1.80 m, weight: 76 kg) were scanned and digitized by a Minolta VI-900 scanner. A profile was added to the back of the forearm in order to reduce drag arising from the base of the solid (Figure 2). The projected area of hand and forearm are respectively 0.018 m 2 and 0.037 m 2 .

Figure 2 Kinematic of the hand-forearm stroke Kinematic parameters were imported into the solver via a data file, which the time step is equal to the physical time-step (2.5.10 -3 s). The path analysis began when the tip of the hand entered into the water and finished when it exited. Kinematic data were measured in the context of swimming. A Vicon ® optoelectronic motion analysis system was used [START_REF] Monnet | Measurement of three dimensional hand kinematics during swimming with a motion capture system: a feasibility study[END_REF]. Eight cameras were placed around a specific pool: two under a bottom window and six in front of a lateral window. The light emitted by the LEDs was reflected back to the sensor by passive markers (fixed on the hand and forearm of the swimmer) covered with specific retroreflective tape. The 3D coordinates of the markers were reconstructed based on the 2D image of two or more cameras whose locations and internal parameters are known. Five reflective markers were fixed on the segment (Figure 2): the finger tip (FT), second (M2) and fifth (M5) metacarpophalangeal joints, radial styloid (RS) and ulnar styloid (US). Velocity (of the FT point) and hand orientation were expressed in a fixed pool-centric reference system: (X0, Y0, Z0) the X0-axis points in the swimming direction, the Y0-axis is perpendicular to the swimming direction pointing from the right to the left of the subject, and the Z0-axis points vertically upward. A local system of reference (Xh, Yh, Zh), linked to the hand-forearm, allowed to calculate the hand orientation (Figure 3). Xh is the axis which passed through FT and the middle (RS, US). Zh was perpendicular to the (Xh, M2M5) plane and Yh was perpendicular to Xh and Zh. The orientations of the hand were respectively, the flexion, pronation-supination, and abduction-adduction of the hand-forearm (Figure 3).

Hand-forearm kinematic parameters

According to [START_REF] Samson | Kinematic hand parameters in front crawl at different paces of swimming[END_REF], some kinematic parameters are directly responsible forces generated by water on the hand and forearm. Five main parameters were defined: velocity, acceleration, and orientation of the hand in the absolute reference system, and the angles of attack and sweepback. In the context of this simulation, these parameters are defined in figure 5.

Figure 5

Analysis methodology of the simulation

The aim of the study is twofold:

1-Study the propulsive forces generated by the hand and forearm during the aquatic stroke and explain it in relation to the kinematic parameters previously defined.

2-Compare the results of the simulation to other similar studies.

For this we will rely on the pressure, force and vorticity field calculated. The propulsive forces result from the action of the flow of water on the surface of the hand and the forearm. Often the resulting forces are expressed in a local reference frame of the flow, making possible a discussion of contributions in lift and drag [START_REF] Schleihauf | A hydrodynamic analysis of swimming propulsion[END_REF][START_REF] Gourgoulis | The influence of the hand's acceleration and the relative contribution of drag and lift forces in front crawl swimming[END_REF]. We will make the choice in this study to decompose the hydrodynamic forces into a pressure component (normal to the wall) and a friction component (tangential to the wall) which represent the two contributions of the surface forces [START_REF] Jr | Fundamentals of aerodynamics[END_REF]. The aim is in particular to discuss the resulting forces directly from the orientations of the hand and the forearm. Thus drag and lift forces and pressure and friction forces do not express in the same reference frames making their comparison difficult.

Results

The main analysis is to jointly study the spatio-temporal evolution of kinematic hand parameters (Figure 5) and the propulsive forces generated by the hand and forearm (Figure 6). The study combines analysis by phase and the comparison with other studies.

Figure 6 Entry and stretch

During the entry and stretch phase (from t = 0 s to t = 0.12 s, Figure 6), the total propulsive force is negative (close to -5 N), and represents the drag force generated by the hand-forearm. This product with low angles of attack, the hand palm directed towards the pool bottom (Figure 4 b,d), which generates few forces.

Downsweep to catch

Propulsion begins during the downsweep to catch phase, before the foremost point of the path (at t = 0.26 s). A first peak force is reached at t = 0.32 s (101 N), with the forearm as the main contribution (Figure 6). At this time, there is a significant deceleration of the hand (for the X0 component, which passes from -23 to 6 m.s -2 between t = 0.24 s and t = 0.28 s, Figure 4 b), and a rapid increase of the angle of attack (from 25° at t = 0.26 s to 50° at t = 0.32 s, Figure 4 c). This increase can be explained by a translational movement of the hand-forearm (due to the stagnation of the angle of flexion close to 60 ° between t = 0.26 s and t = 0.32 s).

In order to better understand the origin of the propulsive forces, a summary of the ressure and friction contributions has been done on the hand and forearm. Results are shown in the Table 1.

Table 1 Results show that the pressure is the main contribution of propulsive forces. This is explained because the flow is turbulent (Re ≈ 10 5 ) and in a separated regime owing to the large angles of attack, the skin friction becomes negligible [START_REF] Dickinson | Unsteady Mechanisms of Force Generation in Aquatic and Aerial Locomotion[END_REF]Anderson 2010, p.79). Thus, the resultant force is approximately perpendicular to the hand-forearm segment. The evolution of the propulsive force during the path can therefore be based mainly on the evolution of the pressure over time on the surface of the hand-forearm segment. A representation of the evolution of these 8 pressures points over time is shown in Figure 7. We chose to subtract the hydrostatic pressure from the static pressure in order to overcome the effect of the water depth to facilitate the analysis of the flow on the surfaces.

Figure 7

It is found that, overall, there are high pressures on the palm of the hand and forearm, and low pressures on the back. It is noted that, at t = 0.26 s, there is an increase of high pressures on the palm of handforearm, and a decrease of low pressures on the back.

To validate the numerical method, these pressure results were compared with experimental study by direct pressure measurement, performed under the same conditions: in crawl and in sprint [START_REF] Toussaint | Pumped-up propulsion" during front crawl swimming[END_REF], but for a different swimmer. In this study, pressure was measured at the shoulder, elbow, wrist, dorsal side of the hand, and palm of the hand using modified disposable blood pressure transducers. Calculation of pressure by numerical simulation was carried out in the same locations: on the elbow, back and palm of the hand (Figure 8). In the figure below, results give pressure difference relative the dorsal side of the hand.

Figure 8 We find that the two curves have two peaks corresponding to the two phases of IN and UP. The maximum values are close, 7000 Pa for our study, and 6500 Pa for the study of [START_REF] Toussaint | Pumped-up propulsion" during front crawl swimming[END_REF].

These two peaks are arranged differently: the largest peak is during the UP phase for [START_REF] Toussaint | Pumped-up propulsion" during front crawl swimming[END_REF], and during the IN phase for our study. These differences may be encountered [START_REF] Maglischo | Swimming Fastest[END_REF]: some swimmers propel themselves further during the IN phase, while others rather during the UP phase.

To go further into the analysis, the results of the propulsive force generated by the hand are compared with two other studies in similar conditions: [START_REF] Kudo | Prediction of propulsive force exerted by the hand in swimming[END_REF] and [START_REF] Sato | Computational Fluid Dynamics Analysis of Hydrodynamic Force Acting on a Swimmer's Hand in a Swimming Competition[END_REF], made respectively from pressure sensors, and numerical simulation from a kinematic measure. The results of both studies and ours are shown in Figure 9. [START_REF] Kudo | Prediction of propulsive force exerted by the hand in swimming[END_REF] computed the propulsive forces exerted by the hand during submaximal swimming pace from a pressure method: a portable data logger with 12 pressure sensors were placed on the hand. Besides, [START_REF] Sato | Computational Fluid Dynamics Analysis of Hydrodynamic Force Acting on a Swimmer's Hand in a Swimming Competition[END_REF] analysed the hydrodynamic forces acting on the hands of the top swimmers at middle-distance pace, by CFD unsteady Navier-Stokes solver based on an unstructured grid method. The hand kinematic parameters were calculated by two video cameras.

Insweep phase

According to the table 1, the IN phase is the most propulsive phase of the aquatic stroke (72.5 N in mean). During this phase, the angle of flexion is between 60° and 110° (Figure 4 b), i.e. the hand is close to the perpendicular of the forward direction. Moreover during this phase, the angles of attack are the highest in the arm's trajectory (between 50° and 70°). The maximum of propulsive force is achieved at t = 0.40 s. This maximum corresponds to a high angle of attack (70°, Figure 4 c), and an angle of flexion equal to 60°. During this phase, the main contribution of forces arises from the hand (57.3 N for the hand and 15.2 N for the forearm, Table 1) which confirms the results of [START_REF] Berger | Hydrodynamic drag and lift forces on human hand/arm models[END_REF]. The end of insweep marks the end of the first peak force. The flow, during this phase, quickly passes from the little finger to the thumb and wrist (sweepback angles of 70 °, 180 ° and 300 °). Therefore, since during this phase the hand-forearm segment is close to the vertical with respect to the forward direction and the pressure contribution is important, the resultant of the hydrodynamic forces is mainly in the forward direction. These results are in agreement with those of [START_REF] Gourgoulis | The influence of the hand's acceleration and the relative contribution of drag and lift forces in front crawl swimming[END_REF] who showed, from drag and lift components that the resultant of the hydrodynamic forces was oriented in the forward direction.

Furthermore, according to this study, it also appears that the resultant forces are perpendicular to the hand of the swimmer, as it has been shown here.

Upsweep phase

A second peak force is located in the UP phase (Figure 6). UP is the second most propulsive phase (30.8

N in mean, Table 1). The beginning of this phase is a transition effected by a change of direction of the hand (at t = 0.54 s, Figure 4 a,c) and a pronation (from 17° at t = 0.54 to 12 ° at t = 0.64 s, Figure 5 c).

In an effort to verify our numerical model further, a comparison with an experimental study [START_REF] Takagi | Unsteady hydrodynamic forces acting on a robotic arm and its flow field: Application to the crawl stroke[END_REF] was made during this IN-UP transition. In this study, measurements were performed for a hand attached to a robotic arm which mimicked a human performing the stroke. The flow fields underwater near the hand were obtained via 2D particle image velocimetry (PIV). Comparison between numerical comparisons made relative to several studies, we can conclude that there is a good agreement with the results of this study.

Limitations of the current study

This new unsteady CFD method opens up many possibilities of studies of the swimming. However, some limitations of this approach must be considered. First, the blocked wrist of the segment does not completely reproduce the kinematics of the forearm. An improvement should be made to provide two degrees of freedom to the wrist. Second, although the eddy-viscosity models gives relevant physical results, it certainly could be improved with second-moment closure models (Reynolds Stress Model, RSM) or Large-Eddy Simulation (LES), or hybrid methods (such as Detached-Eddy Simulation for example). However these models, although constantly improving, are not yet operational for such unsteady flows.

Conclusion

Computational fluid dynamics is a valuable tool in the perspective of a better understanding of propulsive or resistive mechanisms in swimming. The development of computational tools and turbulence models allows researchers to obtain results that can approximate the experimental measurements. Thus the URANS methodology was used to model the flow around the hand and forearm in front crawl swimming. It has been shown that there was a close link between hand kinematic parameters and the propulsive forces generated. This tool thus provides new light on the propulsive mechanisms in swimming and offers some promising developments, especially for optimizing swimming performance from a parametric study.

  inconsistent Reynolds tensor, represents the influence of the background turbulence (modeled) on the filtered field (resolved).

Figure 3

 3 Figure 3

Figure 4 .

 4 Figure 4. Illustration of the hand-forearm segment path with the free surface, as well as the different phases of the aquatic stroke: a) ES phase; b) DC phase; c) IN phase; d) UP phase.

Figure 9

 9 Figure 9

this two studies is based on the vorticity tensor, which is defined as the rotational of the fluid velocity field:

(10) Vorticity (Ω) allows one to know the intensity and the direction of rotation of the vortex structures. The comparison of our results with those of [START_REF] Takagi | Unsteady hydrodynamic forces acting on a robotic arm and its flow field: Application to the crawl stroke[END_REF] is shown in Figure 10.

Figure 10 In Figure 10, we see in both cases a vortex shedded and an induced jet, due to the rapid change of direction. [START_REF] Takagi | Unsteady hydrodynamic forces acting on a robotic arm and its flow field: Application to the crawl stroke[END_REF] shown, from pressure measurements on the surface of the profile, that the jet and vortex shedding was responsible of the increase of propulsive forces at that time.

Between t = 0.62 s and t = 0.68 s, the propulsive forces of hand are high due to the high angles of attack (near 50° at t = 0.6 s, Figure 5 d) and angles of flexion which are close to 110° (Figure 5 c). The decreasing of the propulsive forces from t = 0.68 s (Figure 6) corresponds to the decreasing of the angles of attack which passes from 40° to 15° at this moment (Figure 5 d), as well as the increasing of the angle of flexion which quickly passes from 110° to 140° on the same time (Figure 5 c). So, the palm of the hand is more and more directed towards the surface of the water. At t = 0.72 s, the end of UP phase is the most rearward point of the hand trajectory. The propulsive forces are low since the angle of attack is close to 10 ° (Figure 4 d), and the angle of flexion is close to 140 ° (Figure 4 c). It is the preparation of the release of the hand of the water and the beginning of the recovery.

In summary, there are two force peaks during the IN and UP phases. The hand is the segment that contributes most to the propulsion. The most important propulsive forces are obtained with the higher angles of attack and with angles of flexion close to 90°. Besides, from the three experimental and


Innocent Mufungizi Mihigo

Trésor Hubert Kawaya

Roda Bongeli Mbopanga

Jean Kabulo Banza

Ridi Diakondua Vuilawo

Ruben Loola Loketo

Mapping and Spatiotemporal Dynamics of Land-use and Land-cover Change Based on the USGS and NASA Platforms from MCD12Q1 Landsat Imagery: A Case Study of Kinshasa, DR Congo

Keywords: Broadleaf Forests, Deciduous Broadleaf Forests, Mixed Forests, Closed Shrublands, Woody Savannas, Savannas, Grasslands Chord diagram, Land cover, Land use, MCD12Q1, Model builder. _______________________________________________________________________________________

Kinshasa is a city with very significant urban growth. It is observed that this urban growth is concentrated in a weak area, which requires control and direction of urban expansion. To do this, it is important to understand the dynamics of cover change and land use over time. To do this, we downloaded Landsat images from the MODIS Land Cover Type Product (MCD12Q1) from 2001 to 2022, with a classification of these images based on remote sensing indices that we processed using geographic information system tools. We extracted the following classes: Evergreen

INTRODUCTION

The survival of humanity depends on the natural resources available. It is necessary to improve the management of these resources and strengthen the knowledge on the interactions between humans and the earth through the resources available to them. To understand these interactions, the study of variations of land use/cover therefore becomes a very important source of information because it highlights anthropogenic activities and spatiotemporal change in the environment [START_REF] Sun | Large kernel spectral and spatial attention networks for hyperspectral image classification[END_REF][START_REF] Chen | Temporal and spatial variation of coastline using remote sensing images for Zhoushan archipelago, China[END_REF]. The city of Kinshasa is experiencing very significant urban growth, but urbanization is concentrated in a small area. This situation results from multiple causes which have limited the physical accessibility of the population to land resources: the very strong growth of the population, the saturation of road infrastructure, urban poverty and the low motorization with which it is associated [START_REF] Jean-Pierre | La croissance de l'urbanisation morphologique à Kinshasa entre 1979 et 2015 : Analyse Densimétrique Et De La Fragmentation Du Bâti[END_REF]. DOI: 10.31695/IJASRE.2023.9.11.6

This rapid population growth results in increased demand for social and infrastructure services, education, health and other basic services that are essential to making the city habitable, which poses a challenge questions and concerns. The city of Kinshasa has been officially considered a megacity since 2013, because of the population which was 10 million [START_REF] Ins | Annuaire statistique 2014[END_REF][START_REF] Sosak | Schéma d'orientation stratégique de l'agglomération kinoise et plan particulier d'aménagement de la zone nord de la ville[END_REF]. Unfortunately, the growth observed is often very unorganized and has a direct link with the creation of avenues and informal neighborhoods [START_REF] Watson | The planet city sweeps the poor away,'*: Urban planing and the 21st century urbanisation[END_REF]. These cases are observed in Kinshasa but they are observable in other cities of the Democratic Republic of Congo [START_REF] Useni | Caractérisation de la croissance spatiale urbaine de la ville de Lubumbashi (Haut-Katanga, R.D. Congo) entre 1989 et[END_REF][START_REF] Saint Moulin | Ville et organisation de l'espace en République Démocratique du Congo[END_REF][START_REF] Flouriot | Congo RDC : Population et Aménagement d'un Immense Pays[END_REF]. These problems affect the expanding peripheries through socioeconomic and environmental impacts due to the absence of infrastructure and community services, long journeys to city centers, erosion, flooding, etc. [START_REF] Lelo Nzuzi | Les bidonvilles de Kinshasa[END_REF][START_REF] Angel | Atlas of Urban Expansion (Vol. 1: area and Density)[END_REF][START_REF] Angel | Making room for a planet of cities[END_REF][START_REF] Lelo Nzuzi | Kinshasa : Planification et Aménagement[END_REF].

Land cover is defined by the Food and Agriculture Organization (FAO) as the (bio)-physical coverage of the surface of the emerged land taking into account the type of use or non-use of these lands by humans [START_REF] Jean-Pierre | Occupation du sol. Maison de la télédétection[END_REF]. Land cover includes all physical materials visible on the Earth's surface. The latter can be natural or resulting from anthropogenic activities and take up in their entirety the surfaces occupied by water, trees, grass, paving, bare soil, the surface occupied by buildings, etc. [START_REF] Peter | What is land cover?[END_REF]. Several indices are used, including Normalized difference vegetation index (NDVI) [START_REF] Jia | Mapping Global Distribution of Mangrove Forests at 10-m Resolution[END_REF][START_REF] Zhao | Toward a better understanding of coastal salt marsh mapping: a case from China using dual-temporal images[END_REF], Ratio vegetation index (RVI) [START_REF] Maeda | Can MODIS EVI monitor ecosystem productivity in the Amazon rainforest?[END_REF], Enhanced vegetation index (EVI) [START_REF] Bhatia | Dinoprostone vaginal insert (DVI) versus adjunctive sweeping of membranes and DVI for term induction of labor[END_REF][START_REF] Franch | Remote sensing based yield monitoring: application to winter wheat in United States and Ukraine[END_REF], Difference vegetation index (DVI) [START_REF] Manna | Mapping distribution of Sundarban mangroves using Sentinel-2 data and new spectral metric for detecting their health condition[END_REF][START_REF] Chai | Estimating corn canopy water content from normalized difference water index (NDWI): an optimized NDWI-Based scheme and its feasibility for retrieving corn VWC[END_REF], Normalized difference water index (NDWI) [START_REF] Li | A Genetic Algorithm-based urban cluster automatic threshold method by combining VIIRS DNB, NDVI, and NDBI to monitor urbanization[END_REF][START_REF] Nagy | Wheat yield forecasting for the Tisza River catchment using landsat 8 NDVI and SAVI time series and reported crop statistics[END_REF], Normalized difference building vegetation index (NDBI) [START_REF] Nagy | Wheat yield forecasting for the Tisza River catchment using landsat 8 NDVI and SAVI time series and reported crop statistics[END_REF], Regulated soil vegetation index (SAVI) [START_REF] Zhang | Evaluating ecological health in the middle-lower reaches of the Hanjiang River with cascade reservoirs using the Planktonic index of biotic integrity (P-IBI)[END_REF], Urban building index (IBI-A, IBI-B, IBI) [START_REF] Hou | Marine floating raft aquaculture extraction of hyperspectral remote sensing images based decision tree algorithm[END_REF][START_REF] Alem | Deep learning models performance evaluations for remote sensed image classification[END_REF].

The MODIS Land Cover Type Product (MCD12Q1) provides a suite of science datasets (SDSs) that map global land cover at annual time step for six different land cover legends. The maps were created from classifications of Spectro temporal features derived from data from the Moderate Resolution Imaging Spectroradiometer (MODIS) [START_REF] Mark | User Guide to Collection 6 MODIS Land Cover (MCD12Q1 and MCD12C1) Product[END_REF]. With the aim of studying the dynamics of land cover in Kinshasa (quantitatively and qualitatively evaluating the evolution of land use classes over time) to understand the evolution of the coverage surface of the different classes and their interactions under different methods in the period from 2001 to 2022, we used data acquired by remote sensing that we processed and analyzed using geographic information system techniques and statistical analysis tools in order to be able to give a conclusion of the work. This study contributes to a knowledge base that can help research environmental balance, ecosystem preservation and sustainable development planning by studying interactions between classes.

METHODS

To achieve the objectives set, we downloaded on-site data from United State Geological Survey (USGS) Earth Explorer and National Aeronotics Space Agency (NASA) which presented data from 2001 to 2022 for the area including Kinshasa which is the study area shown in figure 1. For data processing with geographic information system tools, we used the model builder shown here in figure 2. First we use the Extract by mask tool in Spatial Analyst Tools by putting MCD12Q1 as Input raster and the Kinshasa shapefile as Input raster or feature mask data. As input we have a multiband raster file with its attribute values. The mask is specified in the environment parameter when running the Extract by Mask tool, the output raster has cell values only for the area at the intersection of the environment mask and the data in the input mask. The default output format is a geodatabase raster. Secondly, we use the raster obtained from the first step (Output raster) as Input raster in the raster to polygon conversion tool using the values.

The input raster can be any cell size and must be a valid integer raster dataset. We have chosen the attribute field of the input raster dataset which becomes an attribute in the output feature class.

The polygon obtained will finally be used in the dissolve tool using the gridcode to give the area occupied by the land cover classes.

Fig. 2. Model builder of data processing

Data processing was done in ArcGIS using the tools available in ArcToolBox by appying the model builder shown in the figure above. After extraction and processing, we have had the results in figure 3.

DOI: 10.31695/IJASRE.2023.9.11.6 To study the dynamics of land use evolution, it is important to take time and space into account. To do this, we converted the Rasters into polygons based on the values describing the classes. We then used the dissolve tool in ArcGIS based on the grid code to determine the coverage area of each class. The results of this processing phase are presented as percentages by year interval following classes in table 3. To understand the relationships between classes, we used chord diagram (figure 5). A chord diagram is a graphical method of displaying inter-relationships between data radially around a circle. It represents flows or connections between several entities (called nodes) with the relationships between the nodes typically drawn as arcs connecting the data. Here the size of the arc is proportional to the importance of the flow. Chord diagrams are used to visualize data ranging from business use cases to complex scientific data [START_REF] Holten | Hierarchical Edge Bundles: Visualization of Adjacency Relations in Hierarchical Data[END_REF].

This figure shows the relationships between classes. It is observable that certain classes are closely related such as the Woody Savannas and the Savannas. When the first decreases, the second increases and this can be explained by deforestation. It is the same for the Grasslans and Permanent Wetlands, their curves have the same appearance. The regular evolution of Urban and Built-up Lands directly impacts the other classes although on the 21-year scale, it seems small compared to the area studied.

Fig. 1 .

 1 Fig. 1. Location map of the study area

Fig. 3 .Fig. 3 .

 33 Fig. 3. Land cover and land use of Kinshasa (2001-2016)

Fig. 4a .

 4a Fig. 4a. Dynamics of classes taken one by one

Fig. 5 .

 5 Fig. 5. Chord diagram of land cover and land use (%) form 2001 to 2022

Table 1 . Presentation of unprocessed data

 1

	Entity ID	Year	Visualization	Entity ID	Year	Visualization
	MCD12Q1.A20010	-2001		MCD12Q1.A20120	2012
	01.h19v09.061			01.h19v09.061
	MCD12Q1.A20020	2002		MCD12Q1.A20130	2013
	01.h19v09.061			01.h19v09.061
	MCD12Q1.A20030	2003		MCD12Q1.A20140	2014
	01.h19v09.061.2022			01.h19v09.061
	151163011
	MCD12Q1.A20040	2004		MCD12Q1.A20150	2015
	01.h19v09.061			01.h19v09.061
	MCD12Q1.A20050	2005		MCD12Q1.A20160	2016
	01.h19v09.061			01.h19v09.061
	MCD12Q1.A20060	2006		MCD12Q1.A20170	2017
	01.h19v09.061			01.h19v09.061
	MCD12Q1.A20070	2007		MCD12Q1.A20180
	01.h19v09.061			01.h19v09

.061 2018 MCD12Q1.A20080 01.h19v09.061 2008 MCD12Q1.A20190 01.h19v09.061 2019 MCD12Q1.A20090 01.h19v09.061. 2009 MCD12Q1.A20200 01.h19v09.061

  The classes studied are: Evergreen Broadleaf Forests, Deciduous Broadleaf Forests, Mixed Forests, Closed Shrublands, Woody Savannas, Savannas, Grasslands, Permanent Wetlands, Croplands, Urban and Built-up Lands, Cropland/Natural Vegetation Mosaics, Barren and Water Bodies. The description of the classes was made following the description guide presented in table 2.

				2020
	MCD12Q1.A20100	2010	MCD12Q1.A20210	2021
	01.h19v09.061		01.h19v09.061
	MCD12Q1.A20110	2011	MCD12Q1.A20220
	01.h19v09.061		01.h19v09.061

2022

The data thus obtained was processed using geographic information system (GIS) tools, particularly ArcGIS software.

For statistical analyses, we used Excel and OriginPro software. DOI: 10.31695/IJASRE.2023.9.11.6

3. RESULTS AND DISCUSSIONS

Table 2 . MCD12Q1 International Geosphere-Biosphere Programme (IGBP) legend and class descriptions

 2 [START_REF] Mark | User Guide to Collection 6 MODIS Land Cover (MCD12Q1 and MCD12C1) Product[END_REF].

	Name	Value	Description
	Evergreen Needleleaf Forests	1	Dominated by evergreen conifer trees (canopy
			>2m). Tree cover >60%.
	Evergreen Broadleaf Forests	2	Dominated by evergreen broadleaf and palmate
			trees (canopy >2m). Tree cover >60%.
	Deciduous Needleleaf Forests	3	Dominated by deciduous needleleaf (larch) trees
			(canopy >2m). Tree cover >60%.
	Deciduous Broadleaf Forests	4	Dominated by deciduous broadleaf trees
			(canopy
			>2m). Tree cover >60%.
	Mixed Forests	5	Dominated by neither deciduous nor evergreen
			(40-60% of each) tree type (canopy >2m). Tree
			cover >60%.
	Closed Shrublands	6	Dominated by woody perennials (1-2m height)
			>60% cover.
	Open Shrublands	7	Dominated by woody perennials (1-2m height)
			10-60% cover.
	Woody Savannas	8	Tree cover 30-60% (canopy >2m).
	Savannas	9	Tree cover 10-30% (canopy >2m).
	Grasslands	10	Dominated by herbaceous annuals (<2m).
	Permanent Wetlands	11	Permanently inundated lands with 30-60%
			water
			cover and >10% vegetated cover.
	Croplands	12	At least 60% of area is cultivated cropland.
	Urban and Built-up Lands	13	At least 30% impervious surface area including
			building materials, asphalt, and vehicles.
	Cropland/Natural Vegetation	14	Mosaics of small-scale cultivation 40-60% with
	Mo		natural tree, shrub, or herbaceous vegetation.
	saics
	Permanent Snow and Ice	15	At least 60% of area is covered by snow and ice
			for at least 10 months of the year.
	Barren	16	At least 60% of area is non-vegetated barren
			(sand, rock, soil) areas with less than 10%
			vegetation.
	Water Bodies	17	At least 60% of area is covered by permanent
			water bodies.
	Unclassified	255	Has not received a map label because of missing
			inputs.

Table 3 . Percentage of land cover and land use for each class in the study period CLASS

 3

			EVERGREEN BROADLEAF FORESTS	DECIDUOUS BROADLEAF FORESTS	MIXED FORESTS	CLOSED SHRUBLANDS	WOODY SAVANNAS	SAVANNAS	GRASSLANDS	PERMANENT WETLANDS	CROPLANDS	URBAN AND BUILT-UP LANDS	CROPLAND/NATURAL	VEGETATION MOSAICS	BARREN	WATER BODIES
	GRIDC	ODE	2	4	5	6	8	9	10	11	12	13	14		16	17
		2001	1.4835	0.0073	0.0133	0.0053	4.8055	56.2332	33.0503	0.5298	0.0206	3.8292	0.0073		0.0000	0.0147
		2004	1.22035	0	0	0.04587	3.22498	74.140926	17.192015	0.314571	0	3.853343	0		0.0000	0.007945
		2007	0.743076	0.00136	0.002632	0.013443	3.245325	71.760721	20.159899	0.207616	0	3.85997	0		0.0000	0.005959
	PERCENT	2010 2013	0.597613 0.464224	0 0	0.03732 0.036278	0.021154 0.003346	3.445348 3.364137	69.68632 72.393328	22.013064 19.629244	0.31316 0.203675	0.00136 0	3.878703 3.899809	0 0		0.0000 0.0000	0.005959 0.005959
		2016	0.339395	0	0	0.005959	2.386075	68.756595	24.391233	0.212031	0	3.900768	0		0.001986	0.005959
		2019	0.309965	0	0	0.001272	2.398213	68.768172	24.344705	0.241288	0	3.930426	0		0.001986	0.003972
		2022	0.480767	0	0	0	2.091478	67.933466	25.242677	0.259943	0	3.977765	0		0.001986	0.011917
	To understand the evolution of each class, we processed the data in table 3 to produce the following figures:

DOI: 10.31695/IJASRE.2023.9.11.6

  Closed Shrublands: reached a peak in 2004, then a regression until 2022 with a trend of recovery in 2010; -Woody Savannas: We notice a decrease in the coverage area during the study period; -Savannas: this class is the one with the greatest percentage of coverage, we note that it presents a peak in 2004 and the dynamic is regular over time compared to other classes; -Grasslands: This class is second only to Savannas in terms of coverage percentage. It suffered a decline until 2013, then a slight and irregular rise until 2022; -Permanent Wetlands: This class has three trends. It decreases from 2001 to 2007, then shows growth in 2010, then

	decreases in 2013 and finally shows growth until 2022;
	-Croplands: presents a total fall;
	-Urban and Built-up Lands: this class shows regular growth over time;
	-Barren: This class appears in 2016 and presents regular coverage until 2022;
	-Water Bodies: this class decreases until 2007, it remains regular until 2016 to decrease again in 2019 and finally rise
	while remaining regular in 2022.
	Closed Shrublands
	Woody Savannas
	Savannas
	Grasslands
	Permanent Wetlands
	Urban and Built-up Lands
	Water Bodies
	Evergreen Broadleaf Forests

DOI: 10.31695/IJASRE.2023.9.11.6 -

DOI: 10.31695/IJASRE.2023.9.11.6

International Journal of Advances in Scientific Research and Engineering (ijasre), Vol 9 (11), November -2023

Acknowledgments

Not applicable.

Funding

The authors did not receive any funding for this study.

Availability of data and materials

The data used for this study will be made available upon reasonable request.

DOI: 10.31695/IJASRE.2023.9.11.6

CONCLUSION

This study is based on Landsat images from the MODIS Land Cover Type Product (MCD12Q1) with a classification of these images based on remote sensing indices. We found the following land use and cover classes: Evergreen Broadleaf Forests, Deciduous Broadleaf Forests, Mixed Forests, Closed Shrublands, Woody Savannas, Savannas, Grasslands, Permanent Wetlands, Croplands, Urban and Built-up Lands, Cropland /Natural Vegetation Mosaics, Barren and Water Bodies. After calculating the coverage area of these classes, we studied the variation of the percentage of their coverage over time and the relationships existing between these classes. It was noted that the variation in the percentage coverage of certain classes depends on other classes. This is the case of Woody Savannas and the Savannas; Grasslands and Permanent Wetlands. The Urban and Built-up Lands class presents a regular and continuous evolution over time but remains concentrated in one area. Future investigations should focus on data spanning many years with higher resolution. The processing must focus on high-precision algorithms to achieve impeccable results with a minimum of error on cover and land use classes.

SUPPLEMENTARY INFORMATION

Authors' contributions

Declarations

Compliance with Ethical Standards: This article does not contain any studies involving human or animal subjects.

Competing interests: The authors declare that there are no conflicts of interest.


A chain kinematic model to assess the movement of lower-limb including wobbling masses

Keywords: Wobbling masses, Kinematics, Global Optimization, Joint Torques

Soft tissue and bones kinematics can only accurately be distinguished using invasive methods. We propose a non-invasive method to estimate kinematics of both parts in the lower-limb. Modelled as a chain model, including rigid and wobbling bodies, its states were determined using a least-square algorithm. By comparison to a classical model without wobbling bodies, the marker residual was reduced by 20% but the joint kinematics remains unchanged. Average wobbling bodies displacements reached 6.9±3.5° and 6.9±2.4 mm relative to the modeled bones. This original method is a first step to assess the effect of wobbling masses on joint kinetics.

Introduction

In human movement analysis, joint torques are calculated using inverse dynamics in which each body segment is assumed to be rigid. Since body segments can deform [START_REF] Pain | The influence of soft tissue movement on ground reaction forces, joint torques and joint reaction forces in drop landings[END_REF], the rigidity assumption introduces some bias into joint torque estimates. Indeed, two-dimensional computer simulation models with wobbling masses showed that joint efforts (e.g joint torques and joint forces) were 50% lower than those estimated with rigid bodies only [START_REF] Gruber | A comparative study of impact dynamics: wobbling mass model versus rigid body models[END_REF][START_REF] Gittoes | Soft tissue contributions to impact forces simulated using a four-segment wobbling mass model of forefoot-heel landings[END_REF][START_REF] Pain | The influence of soft tissue movement on ground reaction forces, joint torques and joint reaction forces in drop landings[END_REF]. At impact, the deceleration of the wobbling masses was lower and occurred later than the bone deceleration [START_REF] Challis | Soft tissue motion influences skeletal loads during impacts[END_REF]). To our knowledge, computer simulations remain the only method to have assessed the wobbling masses effect on joint dynamics.

Previous models defined each body segment with two rigid parts, namely bone and wobbling masses that move with respect to the bones with three degrees of freedom (DoF), namely translations along the antero-posterior and longitudinal axes and rotation around the medio-lateral axis of the underlying bones [START_REF] Gruber | A comparative study of impact dynamics: wobbling mass model versus rigid body models[END_REF][START_REF] Gittoes | Soft tissue contributions to impact forces simulated using a four-segment wobbling mass model of forefoot-heel landings[END_REF]. However an in-vivo study [START_REF] Gao | Investigation of soft tissue movement during level walking: Translations and rotations of skin markers[END_REF] has shown that skin marker triads mainly rotate around the longitudinal axis of the thigh and shank. Results of simulation models based on planar models could therefore be biased and thus requiring experimental three-dimensional data.

Experimental studies do not separate the bone and wobbling mass kinematics except in intra-cortical pin based experiments. Such experiments have shown that skin markers had a movement composed of a random, marker-specific trajectory (termed as marker deformation), and a common movement called rigid transformation of markers (RTM) caused by the subjacent wobbling masses [START_REF] Leardini | Human movement analysis using stereophotogrammetry -Part 3. Soft tissue artifact assessment and compensation[END_REF]. In-vivo and ex-vivo studies using intra-cortical pins have shown that RTM is greater than marker deformation [START_REF] Andersen | A linear soft tissue artefact model for human movement analysis: proof of concept using in vivo data[END_REF][START_REF] Grimpampi | Quantative assessment of soft tissue artefact propagation to the marker cluster level[END_REF]. Such results suggest that the assessment of the RTM with respect to the bone can be used to determine wobbling mass kinematics.

Experimentally, only intra-cortical pins, in addition to skin markers, could estimate the soft tissue kinematics with respect to the bones. However, they are invasive and provoke pain, especially in tasks involving large impacts with the ground and their use is limited to a few laboratories since their insertion requires sterile conditions [START_REF] Ramsey | Methodological concerns using intra-cortical pins to measure tibiofemoral kinematics[END_REF]. Pins also constrain the wobbling mass motion and sometimes bend due to tendon and muscle contact forces. All these drawbacks have led to the development of mathematical methods to better estimate bone kinematics by minimizing the effect of soft tissue artefacts (STA) [START_REF] Leardini | Human movement analysis using stereophotogrammetry -Part 3. Soft tissue artifact assessment and compensation[END_REF].

Two kinds of mathematical methods have been developed to minimize the effect of STA. Segmental optimization methods consider each body segment independently with six DoF and minimize the deformation of its cluster of markers [START_REF] Alexander | Correcting for deformation in skin-based marker systems[END_REF]Andriacchi 2001, Taylor, Ehrig et al. 2005). Such an assumption can lead to non-physiological joint translation [START_REF] Leardini | Human movement analysis using stereophotogrammetry -Part 3. Soft tissue artifact assessment and compensation[END_REF]. To overcome this limitation, multibody optimization methods include all segments in a chain model with set DoF. The chain model configuration is obtained by minimizing all distances between actual and modelled markers [START_REF] Lu | Bone position estimation from skin marker co-ordinates using global optimisation with joint constraints[END_REF]. The segmental lengths remain constant and "represent the bones" [START_REF] Ojeda | A method to evaluate human skeletal models using marker residuals and global optimization[END_REF]. Some studies have reported better joint kinematics using a multibodyoptimization rather than a segmental optimization approach [START_REF] Klous | Marker-based reconstruction of the kinematics of a chain of segments: a new method that incorporates joint kinematic constraints[END_REF]Klous 2010, Joao, Veloso et al. 2014). Other studies have introduced a soft tissue artefact model in the multibody optimization, which takes into account RTM and does not affect joint kinematics [START_REF] Richard | Introducing a soft tissue artifact model in multi-body optimization: a feasibility study. 3-D Analysis of Human Movement[END_REF]). However, they did not assess the wobbling mass kinematics and their contribution to joint torques remains unknown.

Assuming that a chain model represents the bones, the inclusion of a body representing the wobbling masses in the chain model as in [START_REF] Pain | The influence of soft tissue movement on ground reaction forces, joint torques and joint reaction forces in drop landings[END_REF] should provide a non-invasive method to assess the wobbling mass kinematics. Such a model would be a first step to experimentally investigate the effect of wobbling masses on joint kinetics. The aim of this study was (1) to propose a chain model of the lower-limb including wobbling masses in order to assess their movement. The RTM was considered as the wobbling mass movement and (2) used to determine their main axis of movement. The inclusion of wobbling bodies in a chain model could alter joint kinematics. Hence, the new chain model must estimate the joint kinematics similarly to classical chain models.

Material and Methods

Experimental protocol

Eighteen subjects (23.3±3.4 years old, 1.81±0.11 m, 75.2±10.4 kg, 22.9±2.8 m.kg -2 BMI) took part in this study after giving their informed consent. Eighty-five markers were placed on the lower-limb: pelvis (7), thigh (32), shank (27) and foot (19) (Figure 1a). For the pelvis, the markers were put on the two anterior inferior and superior iliac spines, on the two posterior iliac spines and on the posterior spine of the fifth sacral vertebra. For the thigh, the markers were distributed in six lines of equidistant markers using a tape measure [START_REF] Roosen | Segment-embedded frame definition affects the hip joint centre precision during walking[END_REF] to ensure a reproducible marker placement between subjects. This procedure was repeated for both the shank and the foot.

Please insert Figure 1 here

First, the subjects were instructed to stand still 5 s in an anatomic posture to get the reference position of the markers with respect to their body segment. Subjects then performed functional movements in order to locate the hip and ankle joint centres and the knee flexion axis. Finally, a hopping trial composed of 20 jumps at a self-selected frequency and height of jump was recorded. Subjects received the instruction to bounce with the knees as straight as possible. A force plate was used to detect the landing and take-off times of the jump series. The hopping movement was chosen because it limits flexion at the knee joint and consequently the coupling motions, i.e. translations, adduction and internal rotation [START_REF] Belvedere | Does medio-lateral motion occur in the normal knee? An in-vitro study in passive motion[END_REF]. All marker trajectories were acquired using a 14-cameraoptoelectronic system (Vicon, Oxford, 40 Mpx) sampled at 300 Hz according to the recommendations of [START_REF] Pain | The influence of soft tissue movement on ground reaction forces, joint torques and joint reaction forces in drop landings[END_REF].

Based on all the markers, hip and ankle joint centres were located using the SCoRE algorithm [START_REF] Ehrig | A survey of formal methods for determining the centre of rotation of ball joints[END_REF]. The knee flexion joint axis was determined using the SARA functional method [START_REF] O'brien | Automatic Joint Parameter Estimation from Magnetic Motion Capture Data[END_REF]. The knee joint centre was then defined as the projection on the flexion axis of the mid-point of markers placed on the condyles. The axes of the segmental system of coordinates (SOC) were built to correspond with the anatomical axes (X: medio-lateral, Y: anteroposterior, Z: longitudinal) and their origin was the proximal joint center. The location of each marker is expressed in the SOC of the corresponding body segment and used to determine the state of the chain model.

Lower limb models

Two chain models were built to represent the lower-limb (Figure 1b, Table 1). The first one, called Rigid Model (RM) was a 13 DoF chain model composed of four rigid bodies representing the pelvis, thigh, shank and foot. These bodies were linked by the hip, the knee and the ankle respectively modelled as ball, hinge and ball joints respectively. For the knee, the abduction-adduction and internal-external rotations were considered fixed during the movement. These angles were determined during the static trial using the orientation of the thigh and shank segmental frames for both models.

The second chain model had 25 DoF, named as Wobbling Model (WM). This model added two wobbling bodies to the RM, representing the soft tissues overlying the thigh and the shank with 6 DoF joints for both bodies. The chain models details are summarized in Table 1.

Please insert Table 1 here

The two chain models were reconstructed through a multibody optimization algorithm using a quadratic algorithm similar to [START_REF] Begon | Kinematics estimation of straddled movements on high bar from a limited number of skin markers using a chain model[END_REF] and [START_REF] Fohanno | Improvement of upper extremity kinematics estimation using a subject-specific forearm model implemented in a kinematic chain[END_REF].. The optimal pose (q) of the chain model was obtained for each data frame by minimizing the least squared distance between the actual and the modeled markers in a Newton-Gauss algorithm. The main difference in the model reconstruction between the RM and WM is related to the number of markers involved in the reconstruction process. Indeed, the rigid part of the WM (femur and tibia) is determined using only 36 markers corresponding to those placed on the pelvis, the medial face of the tibia and those placed on the foot. For the WM model, this marker selection means that the marker placed on the thigh and those placed on the posterior face of the shank are used to determine the wobbling mass kinematics. This can be done as the Hessian matrix involved in the Newton Gauss algorithm still remains full ranked.

Data analysis

For the hopping trial, each of twenty jumps was defined as the frame of events between two successive landings. The latter comprised the landing, contact, take-off and mid-aerial phases. Each jump was timenormalized to 200 frames (i.e. 100 frames for contact and aerial phase each). The mean and standard deviation of the joint and wobbling mass kinematics were computed on the determined average-jump.

Model comparisons

The two models were compared using the marker residual to assess the reconstruction quality. This parameter is defined as the mean difference between the measured and the model-determined markers across the entire trial. The difference between the two models was assessed using the Friedman test.

The joint kinematics of RM and WM were compared with regards to joint angles and joint translations.

The RM was used to verify if the joint kinematic obtained was not jeopardized when using the WM model. For each subject, the RMSd was calculated for each joint angle (hip, knee, ankle) and segment translation in the global frame called Tx, Ty and Tz.

A two-way ANOVA with repeated measures was used to compare the models. The second factor corresponded to the four different events of the jump, namely landing, the middle of the contact phase, the take-off and the middle of the aerial phase. If the joint kinematics were similar, then the RTM was analyzed to determine the main axis of the wobbling mass kinematics.

Rigid transformation kinematics

The rotations observed around the medio-lateral, antero-posterior and longitudinal axes were named backward and forward swing, medial and lateral tilt and internal and external twist, respectively. The translations were named as shifts in the antero-posterior, medio-lateral and proximo-distal direction.

The amplitude of rotation and translation on each axis was reported in order to determine the axis with the most motion in rotations and translations. These values were compared using a Kruskal-Wallis test to detect if there were a main axis of translation and rotation. If significant, the Mann-Whitney-Wilcoxon post-hoc test was applied to determine the main axis. All statistics tests were made using IBM SPSS Statistics 19 with α = 0.05, for the kinematic comparisons between the two models and for the determination of the main axis.

RESULTS

The mean RM marker residual is equal to 6.5±2.2 mm whereas the mean WM marker residual is equal to 5.2±0.9mm. The WM model significantly reduced the marker residual of 20% in regards with marker residual obtained using the RM model (p = 0.01).

Joint Kinematics

The RMSd between the two models ranged from 0.4±0.4° to 2.1±1.4° for the joint angles and from 0.7±0.6 mm to 3.5±1.4 mm for the segment translations in the global frame (Table 2). No significant differences were revealed between the two chain models for all joints, rotations and positions involved (p>0.05). For all joints, the mean and standard deviations of joint kinematics for all subjects are illustrated in Appendix A1.

Please insert Table 2 here

Rigid transformation kinematics

The average amplitudes of the thigh wobbling masses rotations and translations among the 18 subjects ranged from 1.2±0.5° to 7.9±3.9° and from 6.2±4.0 to 11.2±4.4 mm, respectively. The statistical analysis demonstrated that the wobbling masses mainly turn around the longitudinal axis and slide along the antero-posterior axis of the thigh (Table 3). Figure 2 shows these results for one subject.

Please insert Table 3 here Please insert Figure 2 here

Regarding the shank, the mean amplitude in rotation of the wobbling masses ranged from 0.9±0.5° to 6.9±3.5° and the translations ranged from 5.4±2.1 mm to 6.9±2.4 mm. The statistical analysis demonstrated that the wobbling masses of the shank mainly turn around the longitudinal axis and slide along the antero-posterior axis (Table 3). These results for another subject are represented in Figure 3.

Please insert Figure 3 here

The mean and standard deviation of the wobbling mass kinematics observed during the average jump among the 18 subjects are available in Appendix A2 and A3 respectively for the thigh and the shank.

Due to high standard deviations, no pattern of the wobbling mass movement among subjects can be defined, suggesting that the movement is subject-specific. Although, for both segments the statistical analysis shows that the wobbling masses mainly turn around the longitudinal axis and mainly slide along the antero-posterior axis.

DISCUSSION

This study is a first attempt to dissociate lower-limb bone and wobbling mass kinematics using skin markers only. This original method consisted in modeling the rigid and wobbling (6 DoF) components of the lower-limb in a multibody chain and estimating the orientation of the femur without taking into account the skin markers located on the thigh. The latter were used to determine the position and orientation of the thigh with respect to the femur only. The main findings are that the WM model gave similar joint kinematics to a classical multibody kinematic approach [START_REF] Lu | Bone position estimation from skin marker co-ordinates using global optimisation with joint constraints[END_REF] and that wobbling masses mainly turn around the longitudinal axis and slide along the antero-posterior axis of both the femur and tibia during hopping. Similarly to the non-invasive methods developed to assess STA (de Rosario, [START_REF] De Rosario | Kinematic description of soft tissue artifacts: quantifying rigid versus deformation components and their relation with bone motion[END_REF][START_REF] Monnet | Assessment of reproducibility of thigh marker ranking during walking and landing tasks[END_REF], the model will have to be validated to estimate its accuracy. Intra-cortical pin based kinematics could be the gold standard measure to validate the model, as only one joint can be study using the bi-planar fluoroscopy [START_REF] Leardini | Human movement analysis using stereophotogrammetry -Part 3. Soft tissue artifact assessment and compensation[END_REF][START_REF] Barre | Soft tissue artifact assessment during treadmill walking in subjects with total knee arthroplasty[END_REF]). However, a restricted sample must be expected [START_REF] Peters | Quantification of soft tissue artifact in lower limb human motion analysis: a systematic review[END_REF].

Model effects

Several studies comparing models demonstrated that the most biofidelic models gave the lowest reconstruction errors (i.e. marker residual) without requiring additional DoF [START_REF] Fohanno | Improvement of upper extremity kinematics estimation using a subject-specific forearm model implemented in a kinematic chain[END_REF][START_REF] Laitenberger | Refinement of the upper limb joint kinematics and dynamics using a subject-specific closed-loop forearm model[END_REF]). In the present study, the WM model also ensured a higher quality of reconstruction by dividing the body segments in two different bodies giving more DoF to the WM model. Its bio-fidelity could be improved using a knee parallel mechanism taking into account the knee joint rotations and translations in function of the flexion-extension joint angle [START_REF] Duprey | Influence of joint constraints on lower limb kinematics estimation from skin markers using global optimization[END_REF]. However the loop closure constraints have to be adjusted so that the final number of DoF at the knee remains equal to 1. The estimation of the femur orientation with respect to the pelvis depends on this assumption, otherwise the Hessian matrix in the Newton-Gauss algorithm is not anymore of full rank.

Model parameters, such as joint type or marker sets associated to each body segment, are known to significantly affect joint kinematics [START_REF] Duprey | Influence of joint constraints on lower limb kinematics estimation from skin markers using global optimization[END_REF]. While no marker was attached directly to the femur, the inclusion of two wobbling bodies representing the soft tissues package had no significant effect on joint kinematics during hopping in a population of 18 young healthy and normal weighted participants. This result fits with the previous work of [START_REF] Richard | Introducing a soft tissue artifact model in multi-body optimization: a feasibility study. 3-D Analysis of Human Movement[END_REF] in which a STA model in a mutlibody optimization was included without seeking to estimate wobbling mass state. Similar to Schulz and Kimmel ( 2010), the WM model also revealed that the tibia and pelvis kinematics are sufficient driving constraints to estimate the femur orientation in dynamics tasks as no significant difference were highlighted in joint kinematics for both studies. In summary, the model proposed in this study can estimate joint kinematics whilst minimizing the reconstruction error in comparison with a classical chain model and thus could be recommended to estimate the wobbling mass movement.

Wobbling mass movement

In hopping with relatively straight legs, wobbling masses mainly rotated around the longitudinal axis and slid along the antero-posterior axis of both the femur and tibia. These results are in agreement with those of [START_REF] Gao | Investigation of soft tissue movement during level walking: Translations and rotations of skin markers[END_REF] who also used skin markers but quantified STA by measuring intermarker translations and rotations during walking trials. This suggests that future computer simulation models built for assessing the wobbling mass effects on joint kinetics should consider these axes of movement which may exert a shear force on the joints. Such aspect will require further evaluation through future research.

While the overall joint kinematics was similar between the two models, a 1-DoF knee joint could bias the wobbling mass kinematics [START_REF] Andersen | Do kinematic models reduce the effects of soft tissue artefacts in skin marker-based motion analysis? An in vivo study of knee kinematics[END_REF]. The values of the rotation and translation of the shank wobbling masses must be considered by the model's assumptions. The angles of abductionadduction and internal and external rotation were modeled as constants and no joint translations allowed whereas these coupling motions vary with the flexion angle [START_REF] Belvedere | Does medio-lateral motion occur in the normal knee? An in-vitro study in passive motion[END_REF]. However the range of flexion-extension (>50°) during hopping stayed far from the values where the knee joint translation is significant according to [START_REF] Belvedere | Does medio-lateral motion occur in the normal knee? An in-vitro study in passive motion[END_REF] and [START_REF] Myers | Measurements of tibiofemoral kinematics during soft and stiff drop landings using biplane fluoroscopy[END_REF].

This reinforces the validity of the model but could limit its application to other movements with greater knee flexion such as running. The tilt and twist movement of the wobbling masses might however be overestimated by this noninvasive approach.

The high between-subjects variability of the wobbling mass movement highlighted that these movements are subject-specific. Furthermore, the difference of mass distribution between subjects (Schinkel-Ivy, Burkhart et al. 2012) and the impact with the ground [START_REF] Pain | Soft Tissue Motion During Impacts: Their Potential Contributions to Energy Dissipation[END_REF] have been considered as the principal sources of variability. A limitation to our study was to not standardize the jump height, which is another source of variability. Additionally, the strategy of muscles recruitment could have affected the wobbling mass kinematics. Accordingly, [START_REF] Boyer | Soft tissue vibrations within one soft tissue compartment[END_REF] observed that subjects with more toned muscles exhibited less soft tissue vibrations. Future studies should correlate wobbling mass kinematics to other biomechanical factors.

The wobbling masses can exhibit different decelerations as a function of their placement on the segment [START_REF] Pain | The influence of soft tissue movement on ground reaction forces, joint torques and joint reaction forces in drop landings[END_REF]. The present study cannot inform which part of the body segment moved the most because the WM chain model associated a single wobbling mass to each bone. The priority for future work would be to introduce the inertial parameters of the different bodies, which remains challenging [START_REF] Challis | Soft tissue motion influences skeletal loads during impacts[END_REF]. Such an estimation could partially rely on the regression equations of [START_REF] Holmes | Predicting in vivo soft tissue masses of the lower extremity using segment anthropometric measures and DXA[END_REF] to differentiate bone and soft tissue masses. Then the model could be applied to assess the contribution of the wobbling masses in the knee shear forces, which are related to ligament injuries, using an inverse dynamic approach.

In conclusion, our original method consists in modeling the rigid and wobbling parts of the segments in a chain model reconstructed using a multibody optimization algorithm. This chain model is able to assess the movement of the wobbling masses in regard to modeled bones without jeopardizing joint kinematics under the assumption of a 1-DoF knee model. The wobbling masses of the thigh and shank mainly turn around the longitudinal axis and slide along the antero-posterior axis of their respective modeled bone. This non-invasive method is a prerequisite to future studies directed at experimentally assessing the effect of the wobbling masses in joint kinetics.

Figure 1

 1 Figure 1. a) Front and back views of the marker placement on the left lower-limb; b) Chain models: Rigid Model (RM) in full line and the two additional wobbling bodies in dotted line which compose the Wobbling Model (WM).

Figure 2 .

 2 Figure 2. Mean kinematic (black line) and standard deviation (grey area) of S5 thigh wobbling masses during the average jump. LA and TO represent landing and takeoff events respectively.

Figure 3 .

 3 Figure 3. Mean kinematic (black line) and standard deviation (grey area) of the S10 shank wobbling masses during the average jump. LA and TO represent landing and takeoff events respectively.

Table 2 .

 2 Mean and standard deviation of RMS differences between the two chain models in terms of joint kinematic and segment positions in the global frames

			Rotations (°)		Translations (mm)
		Flexion-	Abduction-	Internal-	Medio-	Antero-	Vertical
	Joint	Extension	Adduction	External	Lateral	Posterior
	Hip	0.9±0.7	1.0±0.5	2.1±1.4	3.0±1.1	2.4±2.0	0.8±0.7
	Knee	0.9±0.4			1.8±1.9	3.5±1.4	0.7±0.6
	Ankle	0.4±0.4	1.5±1.2	1.2±0.9	2.3±1.4	1.2±0.8	1.6±0.6

Table 3 :

 3 Mean and standard deviation, minimal and maximal amplitudes of the wobbling mass movements for the thigh and the shank wobbling masses. The asterisk highlights the main axis of movement

				Rotation			Translation
						Medio-	Frontward-	Upward-
			Swing (°)	Tilt (°)	Twist (°)	Lateral	Backward	Downward
						(mm)	(mm)	(mm)
	Thigh	Mean Min	1.9±0.8 1.1	1.2±0.5 0.5	7.9±3.9* 4.5	6.2±4.0 2.7	11.2±4.4* 7.4	6.3±3.6 3.6
		Max	3.4	1.8	21.1	11.1	22.7	12.9
	Shank	Mean Min	1.6±0.4 1.4	0.9±0.5 0.4	6.9±3.5* 2.7	5.4±2.1 2.1	6.9±2.4* 5.2	6.1±2.4 3.4
		Max	2.5	1.8	12.4	12.9	11.4	7.6

Acknowledgements

This study was supported by the National Sciences and Engineering Research Council of Canada. The first author received a scholarship from the Foundation of University of Poitiers.

Tables

 Table 1.Characteristics of the RM and WM chain models with body segments and their corresponding markers (see Fig. 1

Appendix

The mean joint kinematics and standard deviation among the subjects obtained with the two models are available in the Figure A1. Figure A1 shows that the two joint kinematics overlap during the average jump. The constant angles for the knee abd-adduction and internal-external rotation are the consequence of the only degree of freedom.

Please insert Figure A1 here Figure A2 and A3 show respectively the thigh and shank wobbling mass kinematics in regard with their respective modeled "bones". In order to illustrate the main axis of the wobbling mass movements in rotation and translation, the same scale was applied for the rotation and translation.

Please insert Figures A2 and A3 here