Chantal Oberson
email: chantal.oberson-ausoni@upmc.fr

Chantal Oberson Ausoni

Pascal Frey
email: pascal.frey@upmc.fr

Geometric algebra for vector field analysis and visualization: mathematical settings, overview and applications

published or not. The documents may come

plex data. Unfortunately, numerical algorithms are too often inadequate to perceive the mathematical properties or the general structure of the objects considered. In this chapter, we will explain how the formal language of geometric algebras may be one of these analysis tools, as it provides a unified framework bringing us closer by the topological aspects of geometrical problems, in a wide range of applications, including scientific visualization. The main strength of geometric algebra lies in the elegance and the generality (ubiquity) of its formulations, which can be injected within the classical Euclidean framework as well as in differential geometry. In this perspective, concepts and ideas introduced should not replace existing theories and tools, but complement them and shed new light on them.

Based on the work of Grassmann, Clifford's geometric algebras, born in the mid 19th-century, consider algebraic operators along with three main products to describe the spatial relations characteristic to geometric primitives in a coordinatefree approach. The many possibilities offered by Clifford algebras and geometric algebras (hereafter denoted GA), and especially their geometrically intuitive aspects, have been emphasized by numerous scientists. For instance, the physicist D. Hestenes has acknowledged their importance to relativistic physics [START_REF] Hestenes | Clifford Algebra to Geometric Calculus: A Unified Language for Mathematics and Physics[END_REF]. Likewise, the mathematicians G.-C. Rota [START_REF] Grosshans | Invariant theory and superalgebras[END_REF], I.R. Porteous [START_REF] Porteous | Clifford algebras and the classical groups[END_REF] and J. Snygg [START_REF] Snygg | A new approach to differential geometry using Clifford's geometric algebra[END_REF], among others, have largely promoted the geometric compactness and simplicity of GA, hence contributing to broaden the field to further applications in computer graphics and scientific visualization.

The next section will briefly present the main concepts and the basic manipulation rules of Clifford and geometric algebras. Then, the specific case of vector fields defined on d-dimensional spaces or on differential manifolds will be addressed in Section 3. In the last section, we will show how geometric algebra can be efficiently used to understand the algebraic structure of vector fields and implemented.

Fig. 1 Sampling a vector field over a cube (left) and summing the trivectors on the unit sphere (right), to compute an approximation of the index, see 4.2 (reprinted from [START_REF] Dorst | Geometric Algebra for Computer Science (Revised Edition): An Object-Oriented Approach to Geometry[END_REF]). Note that such a trivector is a volume in space and not a triple of vectors.

Clifford and geometric algebras

Leaning on the earlier concepts of Grassmann's exterior algebra and Hamilton's quaternions, Clifford intended his geometric algebra to describe the geometric prop-erties of vectors, planes and eventually higher dimensional objects. Basically, Clifford algebra for R n is the minimal enlargement of R n to an associative algebra with unit capturing the metric, geometric and algebraic properties of Euclidean space [START_REF] Gilbert | Clifford Algebras and Dirac Operators in Harmonic Analysis[END_REF]. In general, geometric algebras are distinguished from Clifford algebras by their restriction to real numbers and their emphasis on geometric interpretation and physical applications. Note. Our intent in this section is to give an elementary and coherent account of the main concepts of Clifford and geometric algebras. The reader who is interested in the theoretical aspects of geometric algebras is referred to the textbooks [START_REF] Gilbert | Clifford Algebras and Dirac Operators in Harmonic Analysis[END_REF][START_REF] Hestenes | Clifford Algebra to Geometric Calculus: A Unified Language for Mathematics and Physics[END_REF][START_REF] Hestenes | New foundations for classical mechanics[END_REF], among others. Computational aspects of geometric algebra and its usability in research or engineering applications are discussed in [START_REF] Dorst | Geometric Algebra for Computer Science (Revised Edition): An Object-Oriented Approach to Geometry[END_REF][START_REF] Hildenbrand | Geometric algebra and its application to computer graphics[END_REF]. We privileged a continuous and straightforward digest, deliberately avoiding the conventional succession of definitions and theorems commonly found in most textbooks. Furthermore, most of the concepts in this section are presented in a general setting. The material in this section is intended to be fairly basic but readers unfamiliar with abstract mathematical concepts should skip the formal definition, as well as the advanced concepts in 2.1.2 and 2.2.1.

Clifford algebra

Clifford algebra can be introduced in many ways; the approach we follow here separates the algebraic structure from the geometric interpretation of the product.

Basic notions and definitions

Formal definition. Let V be a vector space over a field K, and let Q : V → K be a quadratic form on V . A Clifford algebra Cl(V, Q) is an associative algebra over K, with identity element 1, together with a linear map i : V → Cl(V, Q) satisfying, for all v ∈ V , the contraction rule i(v) 2 = Q(v)1, such that the following universal property is fulfilled [START_REF] Mahmoudi | Orthogonal symmetries and clifford algebras[END_REF]:

Given any other associative algebra A over K and any linear map

j : V → A such that, for all v ∈ V , j(v) 2 = Q(v)1 A
, there is a unique algebra homomorphism f : Cl(V, Q) → A, for which the following diagram commutes:

V j # # i / / Cl(V, Q) f A
Note that the existence and the uniqueness (up to unique isomorphism) of a Clifford algebra for every pair (V, Q) can be established by considering a quotient algebra of a tensor algebra.

The product defining the Clifford algebra will be called geometric product and denoted as: u v, for u, v ∈ Cl(V, Q) (with a small space between the factors). One usually considers V as a linear subspace of Cl(V, Q), thus dropping the inclusion in the definition of the Clifford algebra, leading u u = u 2 = Q(u). Consequently, the vector space V is not closed under multiplication as, for example, u u is a scalar and not an element of V . The contraction rule also implies that every v ∈ V has an inverse v -1 = v Q(v) , unless Q is degenerate. Intuitive interpretation of the geometric product. One can classically consider the product of two elements a, b ∈ V as the sum of a symmetric and an antisymmetric part:

a b = 1 2 (a b + b a) <a,b> + 1 2 (a b -b a) a∧b .
-In this setting, the symmetric part < a, b > corresponds to the bilinear form φ associated to the quadratic form Q thanks to the polarization identity:

φ (a, b) = 1 2 (Q(a + b) -Q(a) -Q(b)) = 1 2 ((a + b) (a + b) -a a -b b) =< a, b > , this, of course, as a consequence of the contraction rule v 2 = Q(v). When Q is non-degenerate, it is an inner product.
-The antisymetric part a ∧ b has, if non-zero, to be understood as a new entity, that is neither a scalar nor a vector. For Q non-degenerate, the so defined outer product has a very simple interpretation: a ∧ b, for a, b ∈ V , geometrically represents an oriented plane segment, and can be characterized by an algebraic area (the usual area of the parallelogram with the vectors a and b as sides) and the attitude (angular position) of this plane 1 .

The graded Clifford algebra. Consider again the Clifford algebra Cl(V, Q), V and Q like above. We define G 0 as the inclusion of the scalars K in Cl(V, Q). Given an orthonormal basis {e 1 , e 2 , . . . } of V , let G n be the part of Cl(V, Q) generated from the products

n ∏ j=1 e i j , for 1 ≤ i 1 < • • • < i n . The direct sum ∞ n=0
G n is then the graded Clifford algebra. The elements of G n are called n-vectors, where n is the grade.

Elements can be of "mixed grade", like the product a b of two elements in V , which is a sum of a scalar (grade 0) and a bivector (grade 2). A multivector A can be decomposed as a sum A = ∑ ∞ r=0 A r , where A r = A r is of grade r. Extension of the definition of outer product. The outer product of two multivectors A k (grade k) and B ℓ (grade l) is defined as the grade |k + ℓ|part of the product A k B ℓ , writing A k ∧ B ℓ = A k B ℓ k+ℓ . This product extends by linearity on the whole Clifford algebra. For any n ≤ dim(V ), n-blades are defined recursively as outer products of n vectors

a 1 ∧ • • • ∧ a n = (a 1 ∧ • • • ∧ a n-1 ) ∧ a n .
By convention, 0-blades are scalars. A n-blade is a n-vector, but the converse is not true. More precisely [START_REF] Fontijne | Efficient implementation of geometric algebra[END_REF], a sum of two blades A and B is another blade iff they are of the same grade k and share a common factor of grade k -1 or k.

Advanced concepts

Factorization of blades with the geometric product yields two equivalent forms for a blade: one based on the outer product, the other on the geometric product. Actually, for any arbitrary quadratic form Q, given a k-blade A k , it is possible to find an orthogonal basis {v 1 , • • • , v k } of this blade2 . It implies the double formulation

A k =v 1 ∧ v 2 ∧ • • • ∧ v k = v 1 v 2 • • • v k . For example, if a, b ∈ V , with Q(a) non-zero, we have a ∧ b = a ∧ b -<a,b> Q(a) a = a b -<a,b> Q(a)
a . The meet and join are non-linear operations, corresponding to the blade intersection and union. Suppose we have an orthogonal factorization of two blades A and B, i.e., they are given with their orthogonal factorizations A = A ′ C and B = C B ′ , C being the largest common factor. In this very simple case 3

, M = A ∩ B = C and J = A ∪ B = (A ′ C) ∧ B ′ .
Two important involutions are defined on Cl(V, Q): reversion and grade involution. On a r-blade A = (a 1 ∧ a 2 ∧ • • • ∧ a r ), the reversion A † consists of reversing the order of the constitutive vectors (or, because the outer product is antisymmetric on vectors, changing the sign r(r -1)/2 times); the grade involution A # consists of reversing the sign of every constitutive vector:

A † = a r ∧ a r-1 ∧ • • • ∧ a 1 = (-1) r(r-1)/2 a 1 ∧ a 2 ∧ • • • ∧ a r A # = (-1) r A .
The reversion and grade involution extend by linearity on Cl(V, q

): if A = ∑ ∞ r=0 A r , A † = ∞ ∑ r=0 (-1) r(r-1)/2 A r A # = ∞ ∑ r=0 (-1) r A r .
The even (resp. odd) multivectors are the ones with

A # = A (resp. A # = -A).
Using the reversion and the selection of the scalar part • 0 , let us define a bilinear form on Cl(V, Q). On blades A k and B ℓ , we set:

A k * B ℓ =    A † k B ℓ 0 if k = ℓ = 0 A 0 • B 0 if k = ℓ = 0 0 else .
Extending it linearly to multivectors A and B, we obtain the general formula A * B = A † B 0 . Proof of the equivalence between both formulations can be found in [START_REF] Hestenes | Clifford Algebra to Geometric Calculus: A Unified Language for Mathematics and Physics[END_REF], p.13. On vectors, this bilinear form clearly corresponds to the inner product: a * b =< a, b >. When Q is non-degenerate, it is non-degenerate, and it is sometimes called Clifford scalar product.

Geometric algebras

The case V = R n and Q non-degenerate leads to a series of specific definitions and results. As a matter of fact, we have for example:

-Every non-zero blade A r has an inverse

A † r A * A = A † r A † A 0 . -If, in addition, Q is positive definite, then we can define the modulus of element A as |A| = √ A † * A = A † A 0 , since for an element a 1 • • • a r , (a 1 • • • a r ) † (a 1 • • • a r ) = Q(a r ) • • • Q(a 1 ) ≥ 0.
-In R 3 , the existence of an inverse vector has a very clear interpretation. For a given vector v ∈ R 3 and a given scalar a, the equation < v, w >= a defines the affine plane w 0 + v ⊥ . Likewise, given v and a bivector A, the equation v ∧ w = A defines the affine line w 0 + λ v. In both cases, there is no unique solution. However, in the setting of geometric algebra, the equation v w = A leads to the unique solution w = v -1 A (corresponding to the intersection of a plane < v, w >= A 0 and of a line v ∧ w = A 2 ).

Such a Clifford algebra, in the case V = R n and Q non-degenerate, is called geometric algebra. Let (p, q) be the signature of the quadratic form

Q, i.e., Q diag- onalizes in Q(v) = v 2 1 + • • • + v 2 p -v 2 p+1 -• • • -v 2
p+q (Sylvester's law of inertia). We write R p,q for V and Cl p,q for the associated geometric (Clifford) algebra.

Taking a basis {e 1 , e 2 , . . . , e n } of R n , using the element 1 to span the scalars and all products . . e n is called pseudoscalar and is defined to a scalar multiple, since all n-blades are proportional.

r ∏ j=1 e i j for 1 ≤ i 1 < • • • < i r ≤ n (r ∈ N n ) to

Duality and reciprocal frames

The dual A * of a multivector A is defined as

A * = A I -1
n . The duality operation transforms a r-vector A r into an (nr)-vector A r I -1 n ; in particular, it maps scalars into pseudoscalars. The duality relation states (A ∧ B) * = A⌋B * , where ⌋ denotes the left contraction 4 . The inclusion of an element x in a given subspace A specified by a blade A can be defined in two ways:

the direct way:

x ∈ A ⇐⇒ x ∧ A = 0 -the dual way: x ∈ A ⇐⇒ x⌋A * = 0. Given a basis {b 1 , ..., b n } of R p,q , n = p + q, we can define a reciprocal frame {b 1 , ..., b n }, through the formula b i = (-1) i-1 (b 1 • • • ∧ b i-1 ∧ bi ∧ b i+1 • • • ∧ b n ) I -1
n , where I n = b 1 ∧ ... ∧ b n and the ˇ-sign mentions the element removed from the list.

The two basis are mutually orthogonal: < b i , b j >= δ i j . Since the reciprocal of an orthonormal basis is itself, this definition is needed only in non-Euclidean cases. It is also useful in differential geometry.

A vector of R p,q can be written a = ∑ n i=1 a i e i or a = ∑ n i=1 a i e i with a i =< a, e i > and a i =< a, e i >. If we have a multivector basis {e

α |α ∈ {1, • • • , 2 n }, we can also define a reciprocal frame {e α |α ∈ {1, • • • , 2 n }}.

Versors, rotors, spinors and rotations

One of the main features of GA is its ability to deal with the rotations. Indeed, a unique object R can be used to compute the rotation of any subspace X, writing a conjugation with the geometric product:

R(X) = R X R -1 .
The equation x = a x a -1 gives the reflection of an arbitrary vector x along the a-line (a invertible). Its opposite x = -a x a -1 gives the reflection in the dual hyperplane A = a * . Two consecutive reflections form a simple rotation, which can be written as follows:

x" = -b x ′ b -1 = b a x a -1 b -1 = (b a) x (b a) -1 .
It is a rotation of twice the angle between a and b in the plane containing a and b. The element a b is called a 2-versor. In general, a k-versor is a multivector that can be written as the geometric product of k invertible vectors v = v 1 v 2 ... v k . By the Cartan-Dieudonné Theorem [START_REF] Cartan | Lec ¸ons sur la théorie des spineurs[END_REF][START_REF] Dieudonné | Sur les groupes classiques[END_REF], every isometry of R p,q can be reduced to at most n = p + q reflections in hyperplanes. It means that we can write every orthogonal transformation f with a k-versor U (k ≤ n) and the conjugation:

f (x) = (-1) k U xU -1 .
In all spaces of signatures (n, 0), (0, n), (n -1, 1) or (1, n -1), including the Euclidean spaces, every rotation can be written in exponential form 5 :

R(x) = S x S † with S = e 1 2 (i 1 θ 1 +•••+i m θ m ) , i 1 , • • • , i m orthogonal 2-blades .
Note that a rotation of a non-Euclidean space is defined to be an orthogonal transformation of determinant one continuously connected to identity. The element S given by the exponential form of preceding equation is a rotor, i.e., an even versor S satisfying S S † = 1.

A linear map f: V → V can be extended in a function f : Cl(V, Q) → Cl(V, Q) while preserving the outer product:

f (a 1 ∧ a 2 ∧ • • • ∧ a r ) = f (a 1 ) ∧ f (a 2 ) ∧ • • • ∧ f (a r ) .
It is then called an outermorphism. In particular, the reflection of a blade A k in a dual hyperplane a * is (-1) k a A k a -1 and the rotation of a blade by a rotor is RA k R † according to the previous equations for vectors. 5 Quite naturally, the exponential of a blade A is defined with the usual power series ∑ ∞ k=0 A k k! . The additivity exp(A + B) = exp(A) exp(B) is not true in general. The circular and hyperbolic functions of blades are also defined with power series.

Geometric calculus

Differentiation. We consider a finite-dimensional vector space V with quadratic form Q and a multivector-valued function F : U ⊂ V -→ Cl(V, Q). It comes of no surprise that the directional derivative of F in direction r is simply:

F r (x) = lim s→0 F(x + s r) -F(x) s .
This expression will be most of the time written (r * ∇)F instead of F r , expressing the idea of a scalar product between r and the operator ∇, seen as a vector, as will be clearer below. The linearity in r is straightforward; the sum, the geometric product and the grade are preserved. If we want to differentiate a multivector-valued function F : U ⊂ V -→ Cl(V, Q) directly relative to the variable, we consider a base {e 1 , • • • , e m } of V and the coordinate functions of the vector x in this basis x = ∑ n i=1 x i e i . The directional derivatives along the basis directions,

(e i * ∇) = ∂ ∂ x i
, combine into a total change operator 6 as:

∇ = m ∑ i=1 e i (e i * ∇) meaning ∇ F(x) = m ∑ i=1 e i ∂ F(x) ∂ x i .
Note that we also have to define the differentiation from the right, because of the non-commutativity: for a function

F, F(x) ∇ = ∑ m i=1 (e i * ∇)F(x) e i = ∑ m i=1 ∂ F(x) ∂ x i
e i . Thanks to the geometric product, we can write ∇ as: ∇ F = ∇ ∧ F+ < ∇, F >. In the case of a vector-valued function F, we have the usual definitions of the divergence and curl operators:

curl(F) := ∇ ∧ F = 1 2 (∇ F -F ∇) and div(F) :=< ∇, F >= 1 2 (∇ F + F ∇) .
To write the product rule, accents are necessary to specify on what factor the differentiation acts:

∇(F G) = ∇ F G + ∇F G.
The definition of a differentiation with respect to a multivector, for a function

F : U ⊂ Cl(V, Q) -→ Cl(V, Q), is quite straight- forward,
given a reciprocal frame for the whole space Cl(V, Q).

Integration. Consider again a multivector-valued function F; the line integral is

C F(x) dx = lim n→∞ n ∑ j=1 F j ∆ x j , with F j = 1 2 (F(x i ) + F(x i-1 ))
where the chords ∆ x i = x ix i-1 correspond to a subdivision of the curve C. The measure dx is said to be a directed measure, since it is vector-valued. The product between F(x) and dx is the geometric product. If F is vector-valued,

C F(x) dx = C < F(x), dx > + C F(x) ∧ dx . Similarly, if D ⊂ R 2 is a triangulated planar domain, F k is the average of F over the k-th simplex, D F(x) dx = lim n→∞ n ∑ k=1 F k ∆ x k .
The surface measure of the k-th simplex given by vertices x 0 , x 1 , x 2 is

∆ x k = 1 2 (x 0 ∧ x 1 + x 2 ∧ x 0 + x 1 ∧ x 2 ) .
This integral definition can be generalized to higher dimensions [START_REF] Bromborsky | An introduction to geometric algebra and calculus[END_REF].

The fundamental theorem states:

∂V F dS = V F ∇ dX , for a function F defined over a volume V .

Clifford Convolution and Clifford Fourier Transform

For F and H two multivector-valued functions F, H : R m → Cl p,q , the left-and the right-Clifford Convolution of the functions write respectively:

(H * ℓ F)(x) = R m H(x ′ )F(x -x ′ )|dx ′ | , (H * r F)(x) = R m H(x -x ′ )F(x ′ )|dx ′ | .
The quantity |dx| is used to make the integral grade-preserving since dx is a vector within Clifford algebra. Modifying xx ′ into x + x ′ , we get the left-and right-Clifford correlations [START_REF] Ebling | Clifford convolution and pattern matching on irregular grids[END_REF]. The Clifford convolutions generalize the known convolution of scalar-valued functions.

A vector field F can be smoothed through convolution with a scalar field, for example a Gaussian kernel. In the case of two vector fields, the formula for the geometric product leads to the integration of a scalar function < H(xx ′ ), F(x ′ ) > and a bivector function H(xx ′ ) ∧ F(x ′ ) [START_REF] Reich | Analyzing real vector fields with clifford convolution and clifford-fourier transform[END_REF].

In the case of a multivector-valued function F : R 3 → Cl 3,0 , it is possible to define the Clifford Fourier Transform (CFT) of F and its inverse as follows:

F {F}(u) = R 3 F(x) e -2πI 3 <x,u> |dx| , F -1 {F}(x) = R 3 F(u)e 2πI 3 <x,u> |du| .
The function e -2πI 3 <x,u> = cos(2π < x, u >) + I 3 sin(2π < x, u >) is often called Clifford Fourier kernel.

The convolution theorem is also valid for the Clifford Fourier Transform and Clifford convolutions as defined here. For example, using the left convolution,

F {H * ℓ F}(u) := F {H}(u) F {F}(u) .
As mentioned before, the reader willing to get a deeper understanding of the mathematical basics about Clifford algebras and geometric algebras is referred to [START_REF] Gilbert | Clifford Algebras and Dirac Operators in Harmonic Analysis[END_REF][START_REF] Hestenes | Clifford Algebra to Geometric Calculus: A Unified Language for Mathematics and Physics[END_REF][START_REF] Hestenes | New foundations for classical mechanics[END_REF]]. In the next section, we will focus on the analysis of vector fields in the context of GA.

Vector fields in geometric algebra

Our main focus in this paper is the analysis of vector fields, more precisely of steady, linear and non-linear vector fields in Euclidean space and on manifolds. One is classically interested in streamlines, critical points with their classification, separatrices, leading to the topological graph of a vector field. We will show how the analysis of vector fields can benefit from the richer context of geometric algebra.

Vector fields on domains of Euclidean space

Classically, vector fields are mappings of the form v : U ⊂ R n → R n , where U is an open set. Using the notions of GA defined in Section 2.2, if Cl 1 n,0 is the set of 1-vectors of Cl n,0 , a map v : R n -→ Cl 1 n,0 is also a vector field. This definition can be easily extended to bivector, trivector, or spinor fields, for example.

The identification of vector fields (satisfying the Lipschitz regularity condition) with ordinary differential equations dx / dt = v(x) can also be transposed from the classical to the GA setting. The Lipschitz continuity condition can be written in this frame, i.e., there exists a scalar constant K ≥ 0 such that

v(x 1 ) -v(x 2 )) n ≤ K x 1 -x 2 n for all x 1 , x 2 ∈ U .
Furthermore, the defined derivation and integration make it possible to state the existence of an unique solution (streamline or integral curve) through a given point, exactly like in the classical frame. In 2D and 3D, drawing the streamlines is a very classical way to represent a vector field. In order to avoid occlusions and empty areas, departure points (seeds) for these curves are to be placed efficiently.

Let us consider a small example. To a given classical vector field, we can associate curvature and torsion scalar fields: the curvature (resp. torsion) in a point is the curvature (resp. torsion) of the unique streamline in this point [START_REF] Weinkauf | Curvature measures of 3d vector fields and their applications[END_REF]. The curvature field associated to a vector field can be used for the seeding, or can be displayed as a further scalar value in the form of isosurfaces or by color coding. In the GA settings, instead of scalar fields, a curvature bivector field and a torsion trivector field can be defined. Visualizing the curvature bivector along a streamline, we get what is called the Frenet ribbon [START_REF] Benger | Using geometric algebra for visualizing integral curves[END_REF], see figure 2 for such a representation of the vector field.

Vector fields on differential manifolds

Now we turn to vector fields on differential manifolds, having in mind to embed the differential geometry formalism into geometric calculus. For a more detailed presentation of this combined approach, see [START_REF] Hestenes | Clifford Algebra to Geometric Calculus: A Unified Language for Mathematics and Physics[END_REF][START_REF] Snygg | A new approach to differential geometry using Clifford's geometric algebra[END_REF].

Fig. 2 Frenet ribbons constructed from a discrete vector field. The colour encodes the torsion (reprinted from [START_REF] Benger | Using geometric algebra for visualizing integral curves[END_REF]).

In differential geometry, if M is a smooth manifold, a tangent vector in p ∈ M is a derivation, i.e., a linear operator D on C ∞ p (the algebra of germs of smooth functions at p) satisfying the Leibniz rule D( f

• g) = f • D(g) + g • D( f ).
The tangent space of M in p is T p M, the set of such derivations. A vector field is a function assigning to every p ∈ M an element of T p M.

In R 3 , in a more intuitive way, we can imagine giving in each point p of a surface S a vector tangent v(p) to the surface. The link between this v and the associated derivation

D v is the derivative D v ( f )(p) = D f (p)(v(p)).
The operator point of view makes it easier to manipulate vector fields and compose them with other operators. Furthermore a discretization can be made without working with coordinates [START_REF] Azencot | An operator approach to tangent vector field processing[END_REF].

To translate this definition into GA, we give the tangent spaces a Clifford algebra structure. Taking a chart (U, φ ) around p ∈ M, the derivations e p i defined by

e p i ( f ) = ∂ ∂ x i ( f • φ -1 )| x=φ (p)
form a basis for T p M. Forming the blades of these basis vectors, we can build a geometric algebra structure on T p M. With a little more abstraction, a vector field can classically be seen as a section of the tangent bundle, a particular vector bundle: Taking T M to be the disjoint union of tangent spaces on M, T M = ⊔ x∈M T x M, and π : T M → M defined by π(v) = x for x ∈ T x M, we can see M as the base space, T M as the total space and π as the projection, these three elements defining a fibre bundle called the tangent bundle. The section is a continuous map s with π • s = id M , meaning s(x) ∈ T x M, hence what we understand as a vector field. The adding of a geometric algebra structure can be done in the general case of a vector bundle on a manifold with some metrics, using a construction very similar as the one made in 2.1.1: quotienting a tensor algebra with a two-sided ideal.

Scalar fields, vector fields, bivector fields, spinor fields on surfaces, for example, are natural extensions of this definition of vector fields (or can be seen as sections of the Clifford tangent bundle, see above), and, as long as M is simply connected, it is also the case for rotation fields r : M → SO(n), since they can be lifted to spinor fields.

Since every differentiable manifold is locally diffeomorphic to an Euclidean space (via the charts), the existence and uniqueness of streamlines is also granted on manifolds, within or outside GA context.

Critical points, orbits and topological graph

The topological graph is an important tool of analysis: it goes one step further than the streamline representation and decomposes the vector field domain into regions of similar behavior. The critical points and closed orbits (with their type, like defined below) and the separatrices (streamlines or surfaces between areas of different behavior) form the topological graph of the vector field, that eventually describes the underlying structure of the field in a more efficient way as a collection of streamlines. Such a graph does not take into account the norm of the vector field [START_REF] Ebling | Topology based flow analysis and superposition effects[END_REF].

The classification of critical points finds its origin in the theory of dynamical systems. For regular critical points, i.e., for critical points with an invertible Jacobian matrix, a linear approximation of the field can be considered. Studying eigenvalues and eigenvectors of the Jacobian matrix makes the classification possible, provided none of the eigenvalues is pure imaginary. The so-called hyperbolic critical points, satisfying this condition, are isolated and are structurally stable: a small local perturbation does not modify the topology. This justifies the use of the linear approximation to describe the field's behavior around this point. In two dimensions for example, the hyperbolic critical points are sources, sinks, saddles and spirals. Unstable critical points are centers. A similar classification can be done for orbits, according to the derivative of the Poincaré map [START_REF] Asimov | Notes on the topology of vector fields and flows[END_REF]. For non-linear critical points, said to be of higher order, the non-invertibility of the first derivative leads one to consider a higher order Taylor expansion. For the isolated ones, the index 7 can help discriminate critical points of different types. Sometimes this proves insufficient, since two critical points with same index can be of different types. The GA formalism provides an elegant alternative for the computation of the index: for example, Fig. 3 Classification of sectors: hyperbolic, parabolic and elliptic sectors (reprinted from [START_REF] Scheuermann | 17 -topological methods for flow visualization[END_REF]). 7 In 2D, the index corresponds the number of turns the field makes around a critical point.

in 3D,

ind(c) = 1 8 π I 3 B(c) v ∧ dv |v| 3 ,
for v the vector field, c the critical point, B(c) an arbitrary small ball around c [START_REF] Mann | Computing singularities of 3d vector fields with geometric algebra[END_REF]. A corresponding discrete computation will be introduced in 4.2. Unlike the index, the ordered list of all different behavior sectors (i.e., elliptic, hyperbolic and parabolic sectors) makes an unambiguous classification possible [START_REF] Scheuermann | 17 -topological methods for flow visualization[END_REF][START_REF] Henle | A combinatorial introduction to topology[END_REF][START_REF] Firby | Surface topology; 2nd ed. Ellis Horwood series in mathematics and its applications[END_REF][START_REF] Theisel | Topological representations of vector fields[END_REF] (see Fig 3). Next, we turn to a more practical view of geometric algebras, as this chapter is also intended for engineers and practitioners. In particular, we will briefly explain how GA can be implemented and the potential advantages of using Clifford algebra when, for example, dealing with rotations in spaces of high dimensions.

Geometric algebra for computer graphics and visualization of vector fields

Nowadays, geometric algebra is mostly recognized as a promising mathematical concept and is beginning to find broader application. Emerging computer architectures (multicore, many-core, parallel) lead us to believe that the language of GA may find a new playground and evolve towards what Hildenbrand calls Geometric Algebra Computing (GAC) [START_REF] Hildenbrand | Conformal geometric algebra[END_REF]. However, GA is not yet a widespread method in engineering applications, mainly because of two reasons, academic and practical [START_REF] Perwass | Geometric algebra with applications in engineering[END_REF]. On the one hand, GA combines many mathematical concepts that were developed separately over the years and are taught as such in curriculum. On the other hand, most engineering applications in three-dimensional space can be dealt using standard vector and matrix algebra tools. The goal of this last section is to introduce how GA can be used advantageously in computer graphics applications and vector field analysis and visualization.

Geometric algebra for computer graphics

Computer graphics is surely the most obvious field of application of GA. In geometrical applications, operations and transformations are applied on primitives that are combined to represent an object (model). Linear geometric transformations are usually represented using matrices, vectors and scalars. But while 3 × 3 matrices encode the 3D-rotations about an axis through the origin, quaternions are better suited instead, because they are easier to interpret. The quaternion representation of a rotation is a nearly minimal parametrization that requires only four scalars. Given a quaternion, one can easily read off the axis and angle of the rotation, it is not the case with the Euler angles representation. The composition of rotations in quaternion form is faster and more stable numerically (the renormalization is more efficient than with matrices). Furthermore, the interpolation in the set of quaternions H (for example to get an animated view of a rotated object) consists in defining a path on S 3 which is mapped to SO [START_REF] Batard | Clifford bundles: A common framework for image, vector field, and orthonormal frame field regularization[END_REF]. The Euler angles parametrization, from the 3-torus to SO(3) is not a global diffeomorphism: the uniqueness breaks at some points (problem known as the gimbal lock). This is why, in graphic libraries such as OpenGL, rotations are given in terms of a rotation axis and a rotation angle and converted internally into rotation matrices.

Note that H forms a subalgebra of a geometric algebra [START_REF] Porteous | Clifford algebras and the classical groups[END_REF]: it is trivially isomorphic to the even algebra Cl + 3,0 (the set of even multivectors of Cl 3,0 ). We can identify the unit and the basis elements i, j, k of H with the unit and the products e 1 e 2 , e 1 e 3 and e 2 e 3 of Cl 3,0 [START_REF] Hestenes | New foundations for classical mechanics[END_REF]. Thus, the aforementioned quaternion representation of a rotation (of angle θ around unit axis u):

q = e 1 2 θ (u x i+u y j+u z k) = cos 1 2 θ + (u x i + u y j + u z k) sin 1 2 θ with R(x) = q xq -1
, can be seen as a rotor of Cl 3,0 . The products and the conjugation on H and on Cl + 3,0 are defined likewise. Similarly, identifying the pseudoscalar I of Cl 2,0 with the imaginary unit i of C, we have an algebra isomorphism between Cl + 2,0 and C. Clearly, geometric algebra exhibits structural simplicity in the formulations, but its naive implementation may be far less efficient than classical analytical geometry implementations, especially for high dimensions. Fortunately, GA expressions can benefit from compilation (e.g. operator overloading) and parallelization techniques (including GPUs) [START_REF] Fontijne | Efficient implementation of geometric algebra[END_REF][START_REF] Hildenbrand | Conformal geometric algebra[END_REF].

In practice, GA has been implemented in two ways. The additive approach encodes each multivector of Cl n,0 with its 2 n coordinates. It leads typically to a O(2 2n ) time complexity for linear operations and for products (inner, outer and geometric), and to a storage complexity in O(2 n ). The multiplicative approach, restricted to blades, stores the coordinates of the unit orthogonal factors in a matrix and the magnitude using a scalar. Although the storage complexity is smaller than in the additive approach, there is still no available strategy for an efficient implementation of addition. Nevertheless, factorization and simplification operations allow a trivial implementation of "meet" and "join" operations.

As suggested by its name, the conformal model Cl 4,1 of R 3 can be used to represent various angle-preserving geometries. In this model, all conformal transformations can be represented by versors, especially the ones preserving the Euclidean distances.

To define the conformal model Cl 4,1 of R 3 , two vectors e + and e -are adjoined to the basis vectors e 1 , e 2 , e 3 of R 3 embedded in R 5 . They are chosen to form an orthogonal basis, with e + 2 = Q(e + ) = 1 and e - 2 = Q(e -) = -1. If we define respectively n 0 = 1 / √ 2 (e -+ e + ) and n ∞ = 1 / √ 2 (e -e + ), the new basis {e 1 , e 2 , e 3 , n 0 , n ∞ } is not orthogonal (< n 0 , n ∞ >= -1), but makes intuitive definitions for the model possible. The representation p of a point p b ∈ R 3 in the conformal model is defined by the following mapping:

p = F(p b ) = p b + n 0 + 1 2 p 2 b n ∞ .
The element n 0 has the same translation role as the origin vector e 0 in the homogeneous model. The vector n ∞ represents the point at infinity and the axis of symmetry of the horosphere, the set of elements defined by this equation. The Euclidean distance between two points p b , q b ∈ R 3 is directly proportional to the squared root of the inner product < F(p b ), F(q b ) > of their representations in the model Cl 4,1 . The horosphere is formed of null vectors, i.e., vectors of zero norm, as consequence of the fact that p 2 =< p, p > is proportional to (p bp b ) 2 = 0.

The spheres, planes, circles and lines of R 3 can be expressed in the conformal model space Cl 4,1 with two different conditions, using the inner or the outer product. For the example of the sphere:

-A sphere S(a b , r) centered in a b , with radius r corresponds to: Since any vector x ∈ Cl 4,1 can be written x = F(a b ) ± 1 / 2 r 2 n ∞ , for an a b ∈ R 3 and a r ∈ R, the building blocks of Cl 1,4 are spheres, points (spheres with radius zero) and imaginary spheres (spheres with imaginary radius). The reflection in an hyperplane corresponds to a conjugation by a vector in Cl 4,1 . To the other transformations, translations, rotations and scalings, correspond rotors in exponential form (e.g. T = e -1 / 2 t b n ∞ for the translation of vector t b ). All orthonormal transformations can be expressed by rotors, since translations enjoy this property.

S(a b , r) = F(a b ) -1 / 2 r 2 n ∞ ∈ Cl 4,

Geometric algebra for the visualization of vector fields

For the sake of clarity, we restrict ourselves here to 2D and 3D vector fields or vector fields defined on surfaces embedded in R 3 . The objective is to show that GA allows one to perform the local analysis of the fields using differential geometry in a rather classical way, but offers more flexibility and efficiency when identifying the global structures.

With vector data defined at the vertices of a simplicial triangulation T h or of a regular sampling (Cartesian grid), discrete equivalents of geometric and topological entities (e.g. curve, ball) are needed, as well as interpolations, giving vector values at arbitrary locations. This can be achieved in several ways but requires special attention to avoid ambiguous or non-conformal situations [START_REF]Finite metric spaces, combinatorics, geometry and algorithms[END_REF].

To compute the topological index in 2D, we recast the formulation given in Section 3.3 in a discrete setting [START_REF] Henle | A combinatorial introduction to topology[END_REF]. Let B(c) denote a closed polygonal curve around the critical point. For every couple of neighbor vertices (p 1 , p 2 ), form the bivector 1 / 2 ( ṽ(p 1 ) ∧ ṽ(p 2 )) with the values of the normalized vector field ṽ = v / v . The sum of all bivectors, divided by the volume of the unit disk π, will give an approximation of the winding number of v on the curve, which is in turn an approximation of the index of v in c.

It can be shown that two closed polygonal curves discretizing the same underlying continuous curve lead to the same winding number, as long as they are ε-dense (i.e., any point of the continuous curve between two neighbors will be within εdistance of both neighbors). In a continuous setting, the index of a critical point is well defined as the winding number of every circle containing this only critical point, since a nonvanishing vector field v in the interior of a closed path γ implies a zero winding number of v on γ.

A similar computation can be done for 3D vector fields, on a triangulated surface around the critical point (see figure 1). For a triangle of neighbor vertices (p 1 , p 2 , p 3 ) on this surface, form the trivector 1 / 6 ( ṽ(p 1 ) ∧ ṽ(p 2 ) ∧ ṽ(p 3 )) with the values of the normalized vector field. The normalization factor is 4 3 π for the volume of the unit ball [START_REF] Dorst | Geometric Algebra for Computer Science (Revised Edition): An Object-Oriented Approach to Geometry[END_REF][START_REF] Mann | Computing singularities of 3d vector fields with geometric algebra[END_REF]. For a vector field on a surface, the computation is less straightforward than in 2D, since vectors should be projected on a plane, before the sum is computed.

A common aforementioned technique in visualization is to integrate the vector field along a curve, the integral line (or streamline in a fluid). Given a Lipschitz continuous vector field v defined on an open subset U ⊂ R m , one defines curves γ(t) on U such that for each t in an interval I, γ ′ (t) = v(γ(t)). Picard's theorem states that there exists a unique C 1 -curve γ x for each point x in U, so that γ x (0) = x, and γ ′ x (t) = v(γ x (t)) for t ∈ (-ε, +ε). These curves partition the set U into equivalent classes.

Numerically, the discretization of streamlines relies on an integration method; Euler or Runge-Kutta methods are the most common schemes to advance a point along the integral curve given its previous location and a time step δt. Any such method requires to interpolate the field vector at a new location x. The interpolation, defined on classical vector fields using barycentric coordinates, can be written exactly the same way for GA vector fields v : R m → Cl m,0 (m = 2, 3). For example, if x is contained in a simplex then the linear interpolate reads: v(x) = ∑ l i=1 λ i v i , where v i (resp. λ i ) denotes the values of v at the simplex vertices (resp. corresponding barycentric coordinates). Note that the interpolation of a vector field v defined on a triangulated surface S is not straightforward, since the interpolated vectors need to be defined in the tangent planes.

Not every characteristic of the field lies in the topological graph: features such as vortices, shear zones, shock waves, attachment lines or surfaces are not captured in this description and are very important elements to specify the structure of a vector field. The computation methods reviewed in [START_REF] Post | The state of the art in flow visualisation: Feature extraction and tracking[END_REF] to extract features in vector fields are presented in the classical frame but can be extended naturally to the GA frame. Several scalar fields deliver information on the presence of vortices: the vorticity magnitude, the helicity magnitude, the pressure for example. For instance, the vorticity is exactly half of the curl defined in GA.

In some specific situations, the vector field may exhibit local patterns with repetitions over the domain. Their localization would help to apprehend the overall structure of the field. For example, in 2D, we could look for the repetition of singularities like monkey saddles, zones with axis drain, or S-shaped zones. The following approach is inspired by image processing.

Correlation. Given a 2D (resp. 3D) pattern, i.e., a vector field defined on a small square (resp. cubic) domain, we can compute the Clifford correlation (introduced in 2.2.4) between this pattern and a vector field. At each point of the domain, this function gives the similarity of the vector field (in the neighborhood of this point) with the given pattern [START_REF] Ebling | Clifford fourier transform on vector fields[END_REF]. The correlation implies a convolution (quadratic complexity), which can be replaced, via Clifford Fourier Transform, by a multiplication (linear complexity) in the frequency domain. Furthermore, since the 3D CFT can be written as a sum of four complex Fourier transforms through the identification of the pseudoscalar I 3 with the imaginary unit i, Fast Fourier Transforms can be used. However, the main drawback of this method is related to the necessity to check the presence of a given pattern in all positions, for many scales and in many orientations, or the search of the pattern will not be complete.

Invariants. Suppose that we have again a particular feature (patch) we want to identify in a given vector field. Let us attribute values to the different patches through a mapping. Such a mapping, if it exhibits rotation, translation and invariance is called shortly RTS-invariant. If it is, for example, not rotation invariant, then its value has to be computed for all rotated variants of the patch of interest. A family of RTS-invariants and non-redundant moments of order ≤ d [START_REF] Schlemmer | Moment invariants for the analysis of 2d flow fields[END_REF] can be built for 2D scalar and vector fields, using the complex numbers to get a nice formulation of the rotation invariance in the equations. On figure 4, showing a 2D swirling jet flow dataset, the occurences of a given S-shaped pattern can be seen, as obtained by this method. For 3D scalar functions, one of the ways of defining such moments is to use the spherical harmonic functions as building bricks. To extend to 3D vector fields, complex numbers are no help anymore, and quaternions generate a dimension 4 algebra. If the nice formulation of rotations in Clifford algebra and the existence of a product of vectors seems to pave the way for this generalization, the defining of building bricks (perhaps with the spherical vectorial harmonics) for the moments is the first difficulty, followed by the formulation of a rotation invariance condition. To our knowledge, the extension has not been written yet.

Several alternatives to moments as RTS-invariants are defined in literature. For example, the harmonic power spectrum and harmonic bispectrum defined in [START_REF] Fehr | Local rotation invariant patch descriptors for 3d vector fields[END_REF] for 3D vector fields rely on spherical vectorial harmonics. The theory is explained in the classical frame, using representation theory, but possibly further invariants could be defined and a substantial of clarity could be achieved if using GA.

Heat equation. On a Riemannian manifold M, consider the Clifford bundle obtained from the tangent bundle. The Riemannian metric g i j (p) =< e p i , e p j >, positive definite, leads to Euclidean tangent spaces. Let us define now a connection on the manifold ∇ E compatible with the metric (for example the Levi-Cevita connection) and extend it as ∇ C to the Clifford space such that it preserves the graduation, we define a generalized Laplacian as follows:

∆ C = ∑ i j g i j (∇ C e i ∇ C e j -∑ k Γ k i j ∇ C e k ) .
Considering the heat equation ∂ s t ∂t + ∆ C s t = 0, with initial condition s 0 = s, associated with these operators, the solution is a regularization of the section s. It can be approximated through the convolution with the heat kernel. Varying the operators (Clifford-Hodge, Clifford-Beltrami), different flows are obtained, leading to different regularizations. This approach was introduced in [START_REF] Batard | Clifford bundles: A common framework for image, vector field, and orthonormal frame field regularization[END_REF], and was applied to reducing noise in color images, see figure 5, but not yet, to the best of our knowledge, as a global approach tool for vector fields. Fig. 5 A color image (left) with the corresponding unit vector field of edge orientations (middle) and a Clifford-Beltrami regularization of this vector field (reprinted from [START_REF] Batard | Clifford bundles: A common framework for image, vector field, and orthonormal frame field regularization[END_REF]).

In addition to regularization, heat kernel signatures, like they are defined for scalar fields [START_REF] Sun | A concise and provably informative multi-scale signature based on heat diffusion[END_REF], could be used to define signatures of vector field patches.

  span the multivectors, the set {1, e 1 , e 2 , . . . , e n , e 1 e 2 , e 1 e 3 , . . . , e 1 e 2 . . . e n } will form a basis for Cl p,q , with 2 n = ∑ n r=0 n r elements. The element I n = e 1 e 2 .

  1 with p b ∈ S(a, r) ⇐⇒< F(p b ), S(a b , r) >= 0 . -The sphere containing the four points a b , b b , c b , b ∈ R 3 corresponds to the element: S = F(a b ) ∧ F(b b ) ∧ F(c b ) ∧ F(d b ) ∈ Cl 4,1 with p b ∈ S ⇐⇒ F(p b ) ∧ S = 0.

Fig. 4

 4 Fig.4Search for an S-shaped pattern in a 2D swirling jet flow dataset. The original pattern is shown in a green circle, whereas the found occurences are shown in red circles, overlapping for different scales. The method used is the comparison of computed values for a family of moment invariants on the dataset, comparing with the tabulated values for the pattern (reprinted from[START_REF] Schlemmer | Moment invariants for the analysis of 2d flow fields[END_REF]).

The geometric interpretation of the decompostion of the geometric product in outer and inner products will be explained again for V = R 3 at the beginning of

2.2.

A general demonstration (also valid for a degenerate Q) is given for example in[START_REF] Doran | Geometric algebra for physicists[END_REF], page 88. In Euclidean spaces, the well-known Gram-Schmidt orthogonalization can be used.

The dualization introduced in 2.2 makes more general equations for M and J possible.

For two blades A and B of grades a and b, the left contraction A⌋B is A B b-a when a ≤ b, it is zero otherwise. When blade A is contained in blade B, it equals the geometric product A B[START_REF] Dorst | Geometric Algebra for Computer Science (Revised Edition): An Object-Oriented Approach to Geometry[END_REF].

This explains the notation F r = (r * ∇)F for the directional differentiation.

Acknowledgements This work undertaken (partially) in the framework of CALSIMLAB is supported by the public grant ANR-11-LABX-0037-01 overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" program (reference : ANR-11-IDEX-0004-02).


J.-P Dedieu

N Besic

G Vasile

J Mathieu

Y Durand
email: yves.durand@meteo.fr

F Gottardi

DRY SNOW ANALYSIS IN ALPINE REGIONS USING RADARSAT-2 FULL POLARIMETRY DATA. COMPARISON WITH IN SITU MEASUREMENTS

Keywords: dry snow, polarimetry, entropy

In this paper we describe the benefits of RADARSAT-2 in the analysis of temporal changes in polarimetric parameters linked to the snow cover evolution during the winter season. The presented study took place over an instrumented area in the region of French Alps. The focus is set on the dry snow depth retrieval, using an original method based on principal component statistical analysis (PCA) of the polarimetric parameters values. The results obtained by this mean are compared with the network of simultaneous snow in situ measurements. The most thought-provoking result is the strong inverse correlation between the snow depth above the ice crust and the entropy, reflected through the very high coefficient of determination R 2 = 0.8439. In order to justify this observation, we propose an appropriate physical hypothesis.

INTRODUCTION

Snow cover comprises a mixture of ice crystals, liquid water, and air. At temperatures below 0 • , the proportion of liquid water is negligible and this kind of snow is called dry snow. At C band, estimating the location of dry snow covers and the physical properties of the snow pack remains a challenge since dry snow is a low attenuation medium and only affects the amplitude of the backscattered signal to a certain degree [START_REF] Fung | Modeling of microwave emission and scattering from snow and soil[END_REF]. Due to the complexity of scattering, characterized by a wide variety of underlying media, single polarization information at C band is often not sufficient to distinguish between the physical properties of the dry snow volume and the ground. This is particularly true in mountainous environments with strong topographic effects.

A milestone for snow retrieval at C band was reached in Alpine regions with the use of the full polarimetric acquisition mode. Prompted by significant differences in backscattering between the dry and the wet snow [START_REF] Fawwaz | Snowcover influence on backscattering from terrain[END_REF], a new approach for dry snow mapping was developed by Martini et al. [START_REF] Martini | Dry snow discrimination in alpine areas from multifrequency and multi-temporal sar data[END_REF]. They managed to efficiently demonstrate the capabilities of full polarimetry information, provided by the 1994 SIR-C shuttle mission, for observing dry snow covers at L and C bands, under strong topographic conditions. Operational since April 2008, the RADARSAT-2 satellite (5.4 Ghz, C band) now represents a significant and attractive technical advancement in SAR application for snow monitoring thanks to the quadpolarization acquisition mode (HH, HV, VH, VV), shallow incidence angle viewing capabilities and fine resolution mode. Both of latter two appear to be suitable for Alpine environments.

The principal topic of this paper concerns the advantages of RADARSAT-2 in analysis of the snow cover evolution during the winter season. Using the POLSAR data acquired over an instrumented area in the region of French Alps, we analyse the derived polarimetric parameters with respect to the network of simultaneous in situ snow pack measurements. The principal particularity of the presented analysis is the postprocessing of original polarimetric parameters, based on the Principal Component Analysis (PCA). Having a focus set on the dry snow depth retrieval, we compare these, derived polarimetric parameters with the in situ measured snow depth. The most significant observation is the strong inverse correlation between the depth of the snow above the first ice layer and the entropy. Motivated by a very high coefficient of determination of the performed linear regression (R 2 = 0.8439), we propose in this article a physical hypothesis which appropriately relates the entropy and the snow depth.

METHODOLOGY

Total backscattering of the snow pack as a distributed target can be decomposed as the sum of the four following components [START_REF] Besic | Dry snow backscattering sensitivity on density change for swe estimation[END_REF] (Fig. 1):

1. Snow pack surface component, At C band, we can observe a dominance of the snow surface component [START_REF] Fung | Modeling of microwave emission and scattering from snow and soil[END_REF], in case of wet snow, and a dominance of the underlying ground component (3) in case of dry snow [START_REF] Besic | Stochastically based wet snow mapping with sar data[END_REF]. The sensitivity of the polarimetric measurements acquired by RADARSAT-2 could help to retrieve the components ( 2) and ( 4), under certain conditions.

Snow volume component,

It has been demonstrated that relevant polarimetric parameters such as mean alpha angle (α) and entropy (H) can be extracted from the eigenvalue-eigenvector decomposition of the coherency matrix T , and efficiently used to characterize the scattering properties of natural media [START_REF] Cloude | A review of target decomposition theorems in radar polarimetry[END_REF][START_REF] Ferro-Famil | Unsupervised classification of multifrequency and fully polarimetric sar images based on the h/a/α wishart classifier[END_REF][START_REF] Lee | Polarimetric Radar Imaging: from basics to applications[END_REF]. The α angle is the main identification parameter for the nature of the mean scattering mechanism: α = 0 • for surface scattering, 45 • for volume scattering and 90 • for double bounce scattering. The entropy H describes the random behaviour of the global scattering phenomenon. This is due to the fact that entropy does not describe a particular scattering mechanism, but rather characterizes the global scattering process. Indeed, in case of a pure target scattering, H = 0. On the other side, if the scattering is characterized by several, equally influential, scattering mechanisms, then H = 1. In our case study, we used the algebraic H/α decomposition, proposed by Cloude and Pottier in [START_REF] Cloude | An entropy based classification scheme for land applications of polarimetric sar[END_REF]. The polarimetric parameters are obtained using the referent PolSARPro toolbox [START_REF] Dedieu | Interpretation of a radarsat-2 fully polarimetric time series for snow cover studies in an alpine context[END_REF].

Concerning the observed entropy sensibility versus dry snow depth, we propose the following hypothesis as a physical justification. At the considered frequency (f = 5.3 GHz), given the unavoidable presence of ice crust inside the snow pack due to the metamorphism of the snow during the winter season, the reflection of the layer of ice under a pack of dry fresh snow will dominate over the underlying ground backscattering. Generally, refreezing, by forming a snow crust, leads to a very strong volume scattering (large grains representing strong scatterers) [START_REF] Strozzi | Backscattering measurements of alpine snowcovers at 5.3 and 35 ghz[END_REF]. However, if the ice crust (a continuous layer of ice) is formed, it should be rather associated with a strong surface scattering. The support for this conclusion was found in the literature dealing with the passive microwave sensing of a snow cover. Namely, the formation of the crust significantly decreases the measured emissivity, by decreasing the influence of the underlying layers, causing the latter to be ignored in the modelling [START_REF] Matzler | Extension of the microwave emission model of layered snowpacks to coarse-grained snow[END_REF].

Additionally, the inferior roughness of the ice crust surface should rather cause specular reflection of the incident electromagnetic waves than their scattering. This fact, along with the insignificance of the snow surface and volume backscattering components in case of the dry snow, indicates the dominance of the ground-volume interaction component.

As it is seen in Fig. 2, bigger depth of a snow layer (d) means bigger propagation distance (p) of the reflected electromagnetic wave, which is to be rescattered by the volume. This implies bigger probability of interaction with snow particles, therefore bigger probability of recovering the most dominant mechanism. The post-processing used in obtaining the analysed observations is based on the Principal Component Analysis (PCA). Namely, we created a covariance matrix including all the polarimetric parameters versus the target sites of the study, for all dates of images registration. The PCA process was set to reduce the signal noise of the composite components inferring on the polarimetric values (grain size, density, multiple layers). Thus, the data are reorganized using new reference axes where the first component of statistical weight appears to be the major explanatory factor between the physical parameters of the snow pack and the polarimetric descriptors. Then, a correlation table with the snowfield ground measurements was created for the statistical correlation results of the snow depth versus polarimetry.

Fig. 2. Ground-volume interaction

APPLICATION CASE AND RESULTS

The study area is located in the northern French Alps (45 • 09 N, 6 • 10 E), close to the town of Grenoble. This high mountainous region is part of a 1200 km 2 watershed equipped for hydro-power generation ( Électricité de France (EDF), Grenoble) and instrumented by Météo France for avalanches and meteorological forecasting. The RADARSAT-2 imaging mode selected was fine quad-polarization (full polarimetry): swath width 25 × 28 km, Beam FQ19 in single look complex (SLC) mode for sampled pixels. Ground range by azimuth resolution was 7.51 × 4.73 m. Descending orbits (5:42 UT) with right looking direction were selected to be phased with the ground measurements at representative test sites (snow stratigraphy), together with data from EDF and Météo France automatic snow station networks. In regard of the high mountain environment, incidence angle of 39.6 • was selected for the beam configuration. A total of 10 winter images were registered from 2009 to 2011 and 5 dates of dominant dry snow conditions were selected for the study, including 15 sites of depth measurements, located in flat areas. All dry snow depth values measured concern only the section of the snow pack between the air surface (no wetness) and the underlying ground, or the first consistent layer of refrozen snow (solid crust).

RESULTS

After the PCA processing, the correlation table indicates a clear and strong relationship between snow depth measured above the first ice crust and the entropy (H). It is an inverse correlation, being in accordance with the previously described, presumed backscattering physics (Fig. 2). The correlation results, obtained on 5 different dates, at 15 points, are very satisfying with a coefficient of determination for a deduced linear regression R 2 = 0.8439. Entropy values are ranking from 0.50 for large depths (160 to 180 cm) expressing a relatively pure scattering, to 0.75 for smaller dry snow packs (20 to 40 cm) reflecting more complex scattering process (Fig. 3). Fig. 4 describes the temporal evolution of the entropy parameter versus dry snow depth for a set of different sites and dates throughout the winter season (Julian days). Here, we can observe the increase of entropy as we approach the end of the winter season. Given that this temporal evolution is as well characterized by a decrease in dry snow depth,

CONCLUSION

In the French Alps, snow cover is one of the most important components in the run-off modelling and the snowmelt accounts for 50% of the annual stream flow (Rhone river). Snow water equivalent (SWE) characterise the snow pack, as a function of density and depth. Under dry snow conditions, the depth estimation of the snow pack retrieved from the full polarimetry RADARSAT-2 configuration, merged with density network measurements, could offer a precious help for the SWE calculation at the large scale. With a big number of interrelated variables retrieved by the polarimetric decomposition, we observe that the PCA based post-processing of the POLSAR data can reduce the dimensionality of the data set and extract the most important polarimetric information, related to the physical properties of the target. In our case, the entropy parameter seems to be an efficient tool linked to the dry snow depth measurements. Namely, we observed an inverse correlation between the derived entropy and the snow depth above the first strong ice layer. This behaviour was appropriately physically justified.

Fig. 1 .

 1 Fig. 1. Snow backscattering mechanism

Fig. 3 .

 3 Fig. 3. Linear regression: entropy versus dry snow depth

Fig. 4 .

 4 Fig. 4. Temporal evolution of the polarimetric parameter and the dry snow depth

  This study was carried out with the support of the SOAR/CCRS program (#1341), Canada, the PNTS-2009/APOA2 Program France, and the ISIS/Spot-Image contract 2013-731. The project is currently ongoing under the SOAR-EI #5135 and PNTS 2014 #023 programs.