Robert Eymard

Thierry Gallouët

Raphaèle Herbin

RT k mixed finite elements for some nonlinear problems

Keywords: mixed finite element, two phase flow, gradient schemes

RTk mixed finite elements for some nonlinear problems.

Introduction

The numerical solution of environmental underground studies often involves models which require the approximation of linear and nonlinear heterogeneous and anisotropic diffusion operators for general piecewise regular coefficients and on general meshes [START_REF] Douglas | Global estimates for mixed methods for second order elliptic equations[END_REF][START_REF] Arbogast | A nonlinear mixed finite element method for a degenerate parabolic equation arising in flow in porous media[END_REF][START_REF] Amaziane | Numerical modeling of the flow and transport of radionuclides in heterogeneous porous media[END_REF][START_REF] Chavent | Mathematical models and finite elements for reservoir simulation: single phase, multiphase, and multicomponent flows through porous media[END_REF]. A wide number of numerical schemes based on several different approaches have been developed in the last fifteen years to this purpose. An illustration of the variety of these approaches may be found in the two benchmarks which were held in 2008 (two-dimensional case) and in 2011 (3D case) [START_REF] Herbin | Benchmark on discretization schemes for anisotropic diffusion problems on general grids for anisotropic heterogeneous diffusion problems[END_REF][START_REF] Eymard | 3d benchmark on discretization schemes for anisotropic diffusion problem on general grids[END_REF]. The family of gradient schemes was introduced to synthesize Email addresses: robert.eymard@u-pem.fr (Robert Eymard), thierry.gallouet@univ-amu.fr (Thierry Gallouët), raphaele.herbin@univ-amu.fr (Raphaèle Herbin) some of these approaches and was proven to converge for a large number of nonlinear problems [START_REF] Eymard | Gradient schemes for twophase flow in heterogeneous porous media and Richards equation[END_REF][START_REF] Droniou | Gradient schemes: a generic framework for the discretisation of linear, nonlinear and nonlocal elliptic and parabolic equations[END_REF]. This family contains several well-known schemes, such as conforming and lumped conforming schemes, the hybrid mixed mimetic family of schemes [START_REF] Droniou | Gradient schemes: a generic framework for the discretisation of linear, nonlinear and nonlocal elliptic and parabolic equations[END_REF], and some discrete duality finite volume schemes [START_REF] Droniou | A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods[END_REF]. The aim of this paper is to show that it also contains the RT k mixed finite element method [START_REF] Raviart | A mixed finite element method for 2nd order elliptic problems, Mathematical aspects of finite element methods[END_REF][START_REF] Roberts | Mixed and hybrid methods[END_REF] in the case of a linear diffusion problem, and then apply this gradient discretization to two phase flow problems. Let us first recall the primal and dual formulations of the mixed finite element method [START_REF] Durán | Mixed finite elements[END_REF] for the following linear anisotropic diffusion problem. Let Ω be an open bounded connected polygonal subset of R d , d ∈ N , let λ and λ ∈ R, such that 0 < λ ≤ λ and let M d (λ, λ) denote the set of d × d symmetric matrices with eigenvalues in (λ, λ).

Assuming that Λ is a measurable function from Ω to M d (λ, λ), and f ∈ L 2 (Ω), let ū ∈ H 1 0 (Ω) be the solution to:

Ω Λ(x)∇ū(x) • ∇v(x)dx = Ω f (x)v(x)dx, ∀v ∈ H 1 0 (Ω). (1.1)
We now consider a regular simplicial mesh T of Ω. We denote by E the set of edges (in 2D) or faces (in 3D) of the mesh, and E int the subset of internal faces, i.e. such that there exists (K, L) ∈ T 2 such that σ is an interface to K and L, which we denote by σ = K|L. We introduce the following spaces

H div (Ω) = {ϕ ∈ L 2 (Ω) d , divϕ ∈ L 2 (Ω)}, V h = {v ∈ (L 2 (Ω)) d ; v| K ∈ RT k (K), ∀K ∈ T }, (1.2)
V div h = V h ∩ H div (Ω), (1.3)
W h = {p ∈ L 2 (Ω) ; p| K ∈ P k (K), ∀K ∈ T }, (1.4)
M 0 h = {µ : σ∈E σ → R, µ| σ ∈ P k (σ), µ| ∂Ω = 0}, (1.5)
where

• RT k (K) is the Raviart-Thomas space of order k defined on K with di-

mension d k+d d + k+d-1 d-1 ,
• P k (K) is the space of polynomials of d variables of degree k; its dimension is k+d d .

• P k (σ) is the space of polynomials of d -1 variables of degree k; its dimen-

sion is k+d-1 d-1 .
Let us now first recall the primal formulation of the mixed finite element method for Problem (1.1). The primal formulation reads

(v, q) ∈ V div h × W h , (1.6a
)

Ω w(x) • Λ -1 (x)v(x)dx - Ω q(x)divw(x)dx = 0, ∀w ∈ V div h , (1.6b) Ω ψ(x)divv(x)dx = Ω ψ(x)f (x)dx, ∀ψ ∈ W h . (1.6c)
We then introduce the Arnold-Brezzi formulation [START_REF] Arnold | Mixed and conforming finite element methods; implementation, postprocessing and error estimates[END_REF][START_REF] Sayas | From Raviart-Thomas to HDG[END_REF].

(v, q, λ) ∈ V h × W h × M 0 h , (1.7a)
K w(x) • Λ -1 (x)v(x)dx - K q(x)divw(x)dx + σ∈E K σ λ(x) w| K (x) • n K,σ dγ(x) = 0, ∀w ∈ V h , (1.7b)
K ψ(x)divv(x)dx = K ψ(x)f (x)dx, ∀ψ ∈ W h , ∀K ∈ T , (1.7c) σ µ(x) v| K (x) • n K,σ dγ(x) + σ µ(x) v| L (x) • n L,σ dγ(x), ∀σ = K|L ∈ E int , ∀µ ∈ M 0 h , (1.7d)
where n K,σ is the unit normal vector to σ outward K. It is shown in e.g. [START_REF] Arnold | Mixed and conforming finite element methods; implementation, postprocessing and error estimates[END_REF] that the problems (1.6) and (1.7) admit a unique solution, and that the solutions (v, q) to (1.6) and (1.7) are identical. Moreover, let us recall the error estimate [8, Theorem 5.3 p. 39] (due to Brezzi): there exists δ, only depending on λ, λ, on the regularity of the mesh and on Ω such that

q -u L 2 (Ω) + v -Λ∇u H div (Ω) ≤ δ( inf ψ∈W h ψ -u L 2 (Ω) + inf w∈V div h w -Λ∇u H div (Ω) ). (1.8)
In the sequel, we recall the definition of a gradient discretization, from which it is possible to build a gradient scheme for a large class of problems. We then give two constructions of gradient discretization issued from the mixed finite element method, one from the primal formulation and one from the dual formulation, which lead to the same schemes in the case of a linear diffusion anisotropic heterogeneous problem. We then apply these results to the discretization of a nonlinear model of two-phase flow in porous media.

Gradient discretizations for diffusion problems

A gradient discretization D for a second order elliptic problem posed on the domain Ω, with homogeneous Dirichlet boundary conditions on the boundary ∂Ω, is defined by D = (X D,0 , Π D , ∇ D ), where:

• the set of discrete unknowns X D,0 is a finite dimensional vector space on R,

• the linear mapping Π D : X D,0 → L 2 (Ω) is the reconstruction of an approximate function from the discrete unknowns (also often called "lifting operator").

• the linear mapping ∇ D : X D,0 → L 2 (Ω) d is the reconstruction of an approximate gradient from the discrete unknowns; we shall call it "discrete gradient operator". It is chosen such that

• D := ∇ D • L 2 (Ω) d is a norm on X D,0 .
Let us now give the fundamental properties that we seek when designing a gradient discretization (or when recognizing a gradient discretization in an existing scheme) in order to be able to prove the convergence of the resulting gradient schemes for the approximation of linear and nonlinear steady or unsteady diffusion problems [START_REF] Droniou | Gradient schemes: a generic framework for the discretisation of linear, nonlinear and nonlocal elliptic and parabolic equations[END_REF][START_REF] Eymard | Gradient schemes for twophase flow in heterogeneous porous media and Richards equation[END_REF][START_REF] Eymard | Gradient schemes for the Stefan problem[END_REF][START_REF] Eymard | Applications of approximate gradient schemes for nonlinear parabolic equations[END_REF][START_REF] Eymard | Small-stencil 3d schemes for diffusive flows in porous media[END_REF]].

• Coercivity. Let C D be the norm of the linear mapping Π D , defined by

C D = max v∈X D,0 \{0} Π D v L 2 (Ω) v D . (2.1)
A sequence (D m ) m∈N of gradient discretizations is said to be coercive if there exists C P ∈ R + such that C Dm ≤ C P for all m ∈ N.

• Consistency. Let S D be defined by: ϕ ∈ H 1 (Ω) → S D (ϕ) ∈ [0, +∞) with

S D (ϕ) = min v∈X D,0 Π D v -ϕ L 2 (Ω) + ∇ D v -∇ϕ L 2 (Ω) d . (2.2)
A sequence (D m ) m∈N of gradient discretizations is said to be consistent if, for all ϕ ∈ H 1 0 (Ω), S Dm (ϕ) tends to 0 as m → ∞.

• Limit-conformity. Let W D : H div (Ω) × X D,0 → [0, +∞) be defined by

∀(ϕ, u) ∈ H div (Ω) × X D,0 , W D (ϕ, u) = Ω (∇ D u(x) • ϕ(x) + Π D u(x)divϕ(x)) dx. (2.3)
Note that for a conforming finite element method, we have W D (ϕ, u) = 0.

A sequence (D m ) m∈N of gradient discretizations is said to be limit-conforming if, for all sequence u m ∈ X Dm,0 such that u m Dm is bounded, and for all ϕ ∈ H div (Ω), W Dm (ϕ, u m ) tends to 0 as m → ∞.

• Compactness. A sequence (D m ) m∈N of gradient discretizations is said to be compact if, for all sequence u m ∈ X Dm,0 such that u m Dm is bounded, the sequence (Π Dm u m ) m∈N is relatively compact in L 2 (Ω).

Let us illustrate these concepts on Problem (1.1). Let D = (X D,0 , Π D , ∇ D ) be a gradient discretization; then the related gradient scheme for the discretization of (1.1) is to look for u ∈ X D,0 such that

Ω Λ(x)∇ D u(x) • ∇ D v(x)dx = Ω f (x)Π D v(x)dx, ∀v ∈ X D,0 . (2.4)
The coercivity, consistency and limit-conformity properties for a family of gradient discretizations are sufficient to ensure the convergence of Π D u to ū in L 2 (Ω) and that of ∇ D u to ∇ū in L 2 (Ω) d [START_REF] Eymard | Small-stencil 3d schemes for diffusive flows in porous media[END_REF]. The compactness property is only needed for the convergence of gradient schemes in the case of nonlinear problems [START_REF] Droniou | Gradient schemes: a generic framework for the discretisation of linear, nonlinear and nonlocal elliptic and parabolic equations[END_REF][START_REF] Eymard | Gradient schemes for twophase flow in heterogeneous porous media and Richards equation[END_REF][START_REF] Eymard | Gradient schemes for the Stefan problem[END_REF][START_REF] Eymard | Applications of approximate gradient schemes for nonlinear parabolic equations[END_REF].

Primal mixed finite element and gradient discretizations

Here we construct a gradient discretization (in the sense of Section 2) inspired from the primal mixed finite element formulation (1.6) of Problem (1.1). Let W h be defined by (1.4) and let (χ i ) i∈I be a a family of piecewise polynomial basis functions of degree k on each cell of the mesh, spanning W h . Let us define the gradient discretization D = (X D,0 , Π D , ∇ D ) by:

X D,0 = R I , (3.1a)
Π D u = i∈I u i χ i , ∀u ∈ X D,0 , (3.1b) ∀u ∈ X D,0 , Λ∇ D u ∈ V div h , (3.1c
)

Ω v(x) • ∇ D u(x)dx + Ω Π D u(x)divv(x)dx = 0, ∀v ∈ V div h . (3.1d) Remark 3.1.
As we mentioned in the preceding section, a gradient discretization only relies on the definition of discrete operators and therefore, should be problem independent. However we use here the tensor Λ to define the discrete gradient, which enables to ensure the H div conformity in the sense that (3.1c) holds, so that the discrete gradient is problem dependent. Another natural way to obtain a gradient scheme is to construct a gradient discretization without Λ (which is thus problem independent); in this case, we can still prove the convergence, but the approximation is no longer H div conforming; moreover, the resulting approximate gradient is H div consistent, but not Λ∇ D u, and one can expect lower convergence rate in the case of highly anisotropic and heterogeneous problems.

In order for (3.1) to define a gradient discretization, the system (3.1d) should define one and only one ∇ D u, and

• D := ∇ D • L 2 (Ω) d has to be a norm on X D,0 .
The existence and uniqueness of ∇ D u results from the fact that (3.1d) is a square linear system, whose solution vanishes if the right-hand-side vanishes.

The fact that it defines a norm results from the coercivity property shown in the next theorem 3.2.

Let us now recall some known results on the RT k mixed finite element schemes.

The broken Sobolev space H 1 (T )is the set of functions whose restriction to each simplex K of the mesh belongs to H 1 . First recall that, for (V div h , W h ) defined by (1.3)-(1.4) , there exists [8, Lemma 3.5 page 17] an interpolation operator P k defined by

P k : H T = H div (Ω) ∩ (H 1 (T )) d → V div h , ∀p ∈ W h , ∀v ∈ H T , Ω p(x)div(v -P k v)(x)dx = 0, (3.2)
and there exists α > 0, only depending on the regularity of the mesh [8, Theorem

3.1], such that ∀v ∈ H T , v -P k v L 2 (Ω) d ≤ αh( T ∈T v 2 H 1 (T ) ) 1/2 , (3.3)
where h denotes the size of the mesh T . Let us recall how we deduce from the above properties the standard "inf-sup" condition [START_REF] Roberts | Mixed and hybrid methods[END_REF]: let p ∈ W h , let us prolong p by 0 on a ball B with radius R containing Ω. Then there exists

w ∈ H 1 0 (B) such that ∀q ∈ H 1 0 (B), B ∇w(x) • ∇q(x)dx = B p(x)q(x)dx. (3.4)
Moreover w ∈ H 2 (B) and there exists β only depending on d and R such that

w H 2 (B) ≤ β p L 2 (Ω) . (3.5)
Therefore, since ∇w ∈ H T , we have from (3.3)

∇w -P k ∇w L 2 (Ω) d ≤ αhβ p L 2 (Ω) ,
which shows that

P k ∇w L 2 (Ω) d ≤ (2Rα + 1)β p L 2 (Ω) , (3.6)
which leads to the inf-sup condition.

Proof

• Coercivity.

Let u ∈ X Dm,0 (and therefore that p = Π Dm u ∈ (W h ) m ).

Using (3.6), let w ∈ H 1 0 (B) be defined by (3.4), and let v = P k ∇w ∈ V div h . Thanks to (3.2), we get that

p 2 L 2 (Ω) = - Ω p(x)divv(x)dx.
From (3.1), we get

Π Dm u 2 L 2 (Ω) = Ω v(x) • ∇ Dm u(x)dx.
Thanks to (3.6), we then get

Π Dm u L 2 (Ω) ≤ (αD + 1)β ∇ Dm u L 2 (Ω) d ,
which proves the coercivity property.

• Consistency. Let us check the consistency property on the set R = {ϕ ∈

H 1 0 (Ω); there exists f ∈ C ∞ c (Ω) such that ϕ is solution to (1.1)}. Let ϕ ∈ R. Considering Problem (1.6
) with f = -div(Λ∇ϕ), we define u ∈ X Dm,0 by p = Π Dm u and v = -Λ∇ Dm u. Then, we get from (1.8)

Π Dm u -ϕ L 2 (Ω) + Λ∇ Dm u -Λ∇ϕ H div (Ω) ≤ δ( inf ψ∈(W h )m ψ -ϕ L 2 (Ω) + inf w∈(V div h )m w -Λ∇ϕ H div (Ω) ).
Since the right hand side of the above inequality tends to 0 as m → ∞, we obtain that S Dm (ϕ) tends to 0 as m → ∞. The proof of consistency is then

concluded by density of R in H 1 0 (Ω) (see Lemma 3.3 below). • Limit-conformity. Let (u m ) m∈N such that u m ∈ X Dm,0 and ∇ Dm u m re- mains bounded in L 2 (Ω) d as m → ∞. Let ϕ ∈ H div (Ω), and ϕ m ∈ (V div h ) m be an interpolation of ϕ such that ϕ -ϕ m H div (Ω) tends to 0 as m → ∞. Then W Dm (ϕ, u m ) = Ω (∇ Dm u m (x) • ϕ(x) + Π Dm u m (x)divϕ(x)) dx = Ω ∇ Dm u m (x) • (ϕ(x) -ϕ m (x)) + Π Dm u m (x)(divϕ(x) -divϕ m (x)) dx,
thanks to (3.1). Applying the coercivity inequality, we get that the right term of the preceding inequality tends to 0 as m → ∞, which shows the limit conformity of the sequence.

• Compactness. We consider a sequence (u m ) m∈N such that u m ∈ X Dm,0 and ∇ Dm u m remains bounded in L 2 (Ω) d as m → ∞. Then, thanks to the coercivity property, we first extract a subsequence (samely denoted), such that Π Dm u m weakly converges in L 2 (R d ) to some u ∈ L 2 (R d ) (prolonging by 0 outside Ω). Using the limit-conformity, we get that ∇ Dm u m (prolonging by 0 outside Ω) weakly converges in

L 2 (R d ) d to ∇u, which shows that u ∈ H 1 0 (Ω). Let w m ∈ H 1 0 (B) ∩ H 2 (B) (resp. w ∈ H 1 0 (B) ∩ H 2 (B)
) be defined by (3.4) for p = Π Dm u m (resp. p = u). A classical result is that w m converges in H 1 0 (B) to w (from the weak convergence of the gradient and the convergence of its norm).

Applying (3.5), we get that

∇w m H 1 (Ω) ≤ β Π Dm u m L 2 (Ω) . (3.7)
Letting v = P k ∇w m in (3.1), we get

Ω P k ∇w m (x) • ∇ Dm u m (x)dx + Ω Π Dm u m (x)divP k ∇w m (x)dx = 0,
which provides, thanks to (3.2),

Ω P k ∇w m (x) • ∇ Dm u m (x)dx - Ω (Π Dm u m (x)) 2 dx = 0.
Thanks to (3.3) (using (3.7)) and to the convergence of ∇w m to ∇w in L 2 (Ω) d , we get that P k ∇w m converges in L 2 (Ω) d to ∇w. By strong/weak convergence on the first term, we get

lim m→∞ Ω (Π Dm u m (x)) 2 dx = Ω ∇w(x) • ∇u(x)dx = - Ω div(∇w(x))u(x)dx = Ω u(x) 2 dx.
This shows the convergence of Π Dm u m to u in L 2 (Ω), hence concluding the proof of the compactness of the discretization.

Lemma 3.3 (A density result). Let R = {ϕ ∈ H 1 0 (Ω); there exists f ∈ C ∞ c (Ω) such that ϕ solution of (1.1)}. Then R is dense in H 1 0 (Ω).
Proof The mapping T : H 1 0 (Ω) → H -1 (Ω) defined by u → T (u) = -divΛ∇u is continuous and one-to-one thanks to the Lax-Milgram lemma. Therefore the

inverse mapping T -1 is also continuous. Since C ∞ c (Ω) is dense in H -1 (Ω) and R = T -1 (C ∞ c (Ω))
, the conclusion follows. We then have the following property, which establishes the link between the gradient discretization (3.1) and the primal form of the mixed finite element method.

Theorem 3.4. Using the gradient discretization (3.1), the gradient scheme (2.4) for the approximation of Problem (1.1) is equivalent to the primal formulation (1.6) of the mixed finite element method.

Proof Let u ∈ X D,0 be a solution to (2.4). Let us show that (v = -Λ∇ D u, q = Π D u) is the solution of (1.6). We first observe that (3.1d) ensures (1.6b). Let us now consider ψ ∈ W h , which can therefore be written ψ = Π D v, with v ∈ X D,0 .

We then write, thanks to (2.4),

Ω Λ(x)∇ D u(x) • ∇ D v(x)dx = Ω f (x)ψ(x)dx. Setting v = -Λ(x)∇ D u in (3.1d) with u replaced by v, we have Ω v(x) • ∇ D v(x)dx + Ω Π D v(x)divv(x)dx = 0, which implies Ω Λ(x)∇ D u(x) • ∇ D v(x)dx = Ω ψ(x)divv(x)dx.
This completes the proof of (1.6c). Reciprocally, considering the solution (v, q) to (1.6), since q ∈ W h , there exists a unique u ∈ X D,0 such that q = Π D u.

From (1.6b), we get that v = -Λ∇ D u. Following the same computation lines as the beginning of this proof and letting ψ = Π D v for any v ∈ X D,0 , we get that (1.6c) implies (2.4), again using (3.1d) with u replaced by v.

Mixed hybrid formulation and gradient schemes

Here we construct a gradient discretization (in the sense of Section 2) inspired from the dual mixed finite element formulation (1.7) of Problem (1.1). Let W h be defined by (1.4) and let again (χ i ) i∈I be a a family of piecewise polynomial basis functions of degree k on each cell of the mesh, spanning W h . Let M 0 h be defined by (1.5) and let (ξ j ) j∈J be a family spanning M 0 h . Recall that V h is defined by (1.2). Let us define the gradient discretization D = (X D,0 , Π D , ∇ D ) by:

X D,0 = R I∪J , (4.1a)
Π D u = i∈I u i χ i and Γ D u = j∈J u j ξ j , ∀u ∈ X D,0 , (4.1b)
∀u ∈ X D,0 , Λ∇ D u ∈ V h and:

K w(x) • ∇ D u(x)dx + K Π D u(x)divw(x)dx - σ∈E K σ Γ D u(x) w| K (x) • n K,σ dγ(x) = 0, ∀w ∈ V h . (4.1c)
Remark 4.1. In the case k = 0, the gradient discretization (4.1) has the same degrees of freedom as the SUSHI scheme [START_REF] Eymard | Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes SUSHI: a scheme using stabilization and hybrid interfaces[END_REF], which is also a gradient discretization. Nevertheless, the definitions of the discrete gradients are different.

As in the previous section, in order for (4.1) to define a gradient discretization, the system (4.1c) should define one and only one discrete gradient ∇ D u, and

• D := ∇ D • L 2 (Ω) d has

Proof

Let us first denote by T : X D,0 → X D,0 , T (u) = (u i ) i∈I , where we use the gradient discretization D defined by (3.1). The first remark is that the definition of T implies Π D u = Π D T (u) a.e. in Ω. We observe that, selecting

w = Λ∇ D T (u) ∈ V div h ⊂ V h in (4.1c
) and summing on K ∈ T , all the integrals on σ ∈ E int vanish, and we obtain

Ω Λ∇ D T (u)(x) • ∇ D u(x)dx + Ω Π D u(x)div(Λ∇ D T (u))(x)dx = 0. Now using (3.1d) with v = Λ∇ D T (u), we get Ω Λ∇ D T (u)(x) • ∇ D u(x)dx = Ω Λ∇ D T (u)(x) • ∇ D T (u)(x)dx.
We thus obtain, thanks to the Cauchy-Schwarz inequality,

Ω Λ∇ D T (u)(x) • ∇ D T (u)(x)dx ≤ Ω Λ∇ D u(x) • ∇ D u(x)dx,
which leads to

∇ D T (u) 2 L 2 (Ω) d ≤ λ λ ∇ D u 2 L 2 (Ω) d , ∀u ∈ X D,0 . (4.2)
• Coercivity follows from the coercivity of the gradient discretization D and (4.2).

• Consistency. Thanks to [18, Proposition 3.1 p. 15], for any u ∈ X D,0 , there

exists u ∈ X D,0 such that T ( u) = u and ∇ D u = ∇ D u. The consistency of D m
is then a consequence of that of D m .

• Limit-conformity. Similarly to the study of the limit-conformity of D m , we get

W Dm (ϕ, u m ) = Ω ∇ Dm u m (x) • ϕ(x) + Π Dm u m (x)divϕ(x) dx = Ω ∇ Dm u m (x) • (ϕ(x) -ϕ m (x)) + Π Dm u m (x)(divϕ(x) -divϕ m (x)) dx,
thanks to (4.1c) and since the terms on σ ∈ E int vanish, using that ϕ m ∈ (V div h ) m . So the conclusion is identical.

• Compactness. This property is again an immediate consequence of the compactness of the gradient discretization D and (4.2). Proof Let u ∈ X D,0 be a solution to (2.4). Let us show that (v = -Λ∇ D u, q = Π D u, λ = Γ D u) is the solution of (1.7). We first observe that (4.1c) ensures

(1.7b). Let us now consider ψ ∈ W h and µ ∈ M 0 h . Consider a particular K ∈ T , and take in (2.4) a function test v ∈ X D,0 such that Π D v| K = ψ| K , Π D v| L = 0 for all L ∈ T \{K} and Γ D v = 0. We remark that, thanks to (4.1c), the support of ∇ D u is also reduced to K, and that we may write, thanks to (2.4),

K Λ(x)∇ D u(x) • ∇ D v(x)dx = K f (x)ψ(x)dx.
Setting w = v in (4.1c) where u is replaced by v, we get, using

Γ D v = 0, K v(x) • ∇ D v(x)dx + K Π D v(x)divv(x)dx = 0, which implies K Λ(x)∇ D u(x) • ∇ D v(x)dx = K Π D v(x)divv(x)dx = K ψ(x)divv(x)dx.
This completes the proof of (1.7c). Then, we let v ∈ X D,0 in (2.4) such that Π D v = 0 and Γ D v| σ = λ| σ for a given σ = K|L ∈ E int , and Γ D v| σ = 0 for all σ ∈ E \ {σ}. Again setting w = v in (4.1c) where u is replaced by v, we get,

using Π D v = 0, K v(x) • ∇ D v(x)dx - σ Γ D v(x) v| K (x) • n K,σ dγ(x) = 0, and L v(x) • ∇ D v(x)dx - σ Γ D v(x) v| L (x) • n L,σ dγ(x) = 0.
This implies (1.7d), using (2.4) and the fact that the support of

∇ D v is reduced to K ∪ L.
Conversely, considering the solution (v, q, λ) to (1.7), since q ∈ W h and λ ∈ M 0 h , there exists a unique u ∈ X D,0 such that q = Π D u and λ = Γ D u. From (1.7b), we get that v = -Λ∇ D u. For any v ∈ X D,0 , letting ψ = Π D v and µ = Γ D v, following the same computation lines as the beginning of this proof, we get that (1.7c) and (1.7d) imply (2.4) using (4.1c) where u is replaced by v.

Here again, the fact that we wish to obtain a mixed finite element scheme has led us to a problem dependent discretization. For the linear problem (1.1), this is not a difficulty. However if the diffusion tensor depends on the unknown, it becomes very intricate to ensure the H div conformity, and in fact quite useless since one can get the approximate continuity of the flux from the gradient discretization itself. This line of thought may lead to consider the scheme (4.1),

replacing Λ∇ D u ∈ V h by ∇ D u ∈ V h in (4.1c).

Application to two phase flow in porous media and Richards' equation

We are interested here in the approximation of (u, v), solution to the incompressible two-phase flow problem in the space domain Ω during the time period (0, T ), which was studied in [START_REF] Eymard | Gradient schemes for twophase flow in heterogeneous porous media and Richards equation[END_REF] in the case of a general gradient discretization; here we propose to show that the specific gradient discretization (4.1) obtained from the dual mixed finite element formulation is well adapted to this problem.

The continuous problem reads:

Φ(x)∂ t S(x, p) -div(k 1 (x, S(x, p))Λ(x)∇u) = f 1 , (5.1a)
Φ(x)∂ t (1 -S(x, p)) -div(k 2 (x, S(x, p))Λ(x)∇v) = f 2 , (5.1b)
p = u -v, for ∈ Ω × (0, T ), (5.1c)
considered with the following initial condition: S(x, p(x, 0)) = S(x, p ini (x)), for a.e. x ∈ Ω, (

together with the following non-homogeneous Dirichlet boundary conditions: u(x, t) = ū(x) and v(x, t) = v(x) on ∂Ω × (0, T ).

(5.1e)

In Problem (5.1), u (resp. v) denotes the pressure of the phase 1, called the wetting phase (resp. of the phase 2, which is the non-wetting phase), p is the difference between the two pressures, called the capillary pressure, the saturation of the phase 1 is denoted by S(x, p) (it is called the "water content" in the framework of Richards' equation), and where Φ, Λ, k i , and f i (i = 1, 2) respectively denote the porosity, the absolute permeability, the relative permeabilities, and the source terms, under the following assumptions:

• Ω is an open bounded connected polyhedral subset of R d , d ∈ N and T > 0, (5.2a)
• Φ is a measurable function from Ω to R with

Φ(x) ∈ [Φ min , Φ max ], Φ max ≥ Φ min > 0, (5.2b)
• Λ is a measurable function from Ω to M d , where M d denotes the set of d × d matrices, such that for a.e. x ∈ Ω, Λ(x) is symmetric, and the set of its eigenvalues is included in [λ, λ],

with 0 < λ ≤ λ, (5.2c)
• p ini ∈ L 2 (Ω), (5.2d)
• S(x, q) ∈ [0, 1] for all (x, q) ∈ Ω × R with S(x, q) = S j (q) for a.e. x ∈ Ω j and all q ∈ R, where S j is a non decreasing Lipschitz continuous function with constant L S , (Ω j ) j∈J is a family of disjoint connected polyhedral open sets such that j∈J Ω j = Ω where J is a finite set, (5.2e)

• f i ∈ L 2 (Ω × (0, T )), i = 1, 2, (5.2f) • k i (x, s) ∈ [k min , k max ] for (x, s) ∈ Ω × [0, 1] with k max ≥ k min > 0 and k i (•, s) measurable, k i (x, •) continuous, i = 1, 2, (5.2g)
• ū, v ∈ H 1 (Ω).

(5.2h) Remark 5.1 (Gravity terms). Gravity terms are taken into account in the mathematical analysis performed in [START_REF] Eymard | Gradient schemes for twophase flow in heterogeneous porous media and Richards equation[END_REF]; we choose not to include them here for the sake of brevity and clarity.

Assumptions (5.2) are quite general, except for k min > 0 in Hypothesis (5.2g).

This assumption is needed in the mathematical proof of convergence [START_REF] Eymard | Gradient schemes for twophase flow in heterogeneous porous media and Richards equation[END_REF]. Assumption (5.2e) is compatible with the so-called Van Genuchten model S j (p) = 1/((max(-p, 0)/p j ) nj + 1) mj , with real parameters p j , n j , m j > 0. The hypothesis that the function S(x, p) is defined by given functions in a partition of the domain is classical. Problem (5.1) is considered under the following weak sense.

u -ū, v -v ∈ L 2 (0, T ; H 1 0 (Ω)), p = u -v ; ∀ϕ ∈ C ∞ c (Ω × [0, T )), T 0 Ω -Φ(x)S(
+ Ω Φ(x)S(x, p ini (x))ϕ(x, 0)dx = T 0 Ω f 2 (x, t)ϕ(x, t)dxdt,
where we denote by

C ∞ c (Ω × [0, T )) the set of the restrictions of functions of C ∞ c (Ω × (-∞, T )) to Ω × [0, T ).
Alternately, we also consider the case of Richards' equation, which can be obtained from (5.1) in two ways. One can replace (5.1a) by u(x, t) = ū(x) for (x, t) ∈ Ω × (0, T ), (5.4) or one can replace (5.1b) by v(x, t) = v(x) for (x, t) ∈ Ω × (0, T ).

(5.5)

Then the corresponding weak formulation of the problem is obtained by replacing (5.3a) by (5.4) or (5.3b) by (5.5). Note that this sense is an extension of the condition classically used in hydrogeological studies, which prescribes a constant condition with respect to time and space for the air pressure. As shown in the numerical examples of [START_REF] Eymard | Gradient schemes for twophase flow in heterogeneous porous media and Richards equation[END_REF], this extension allows to use Richards' equation as a good approximation of the full two phase flow problem in other engineering frameworks.

We consider a time interval (0, T ) and (t (n) ) n=0,...,N such that t (0) = 0 < t (1) . . . < t (N ) = T . We then set δt (n+ 1 2 ) = t (n+1) -t (n) , for n = 0, . . . , N -1.

Let us consider the following scheme for the approximation of Problem (5.1) (let us emphasize that this discretization is not based on the global pressure formulation of the problem). We consider the gradient discretization given by (4.1) (we emphasize that we cannot consider the gradient discretization given by (3.1), for an important reason explained below); recall that this discretization is partially problem dependent since it depends on the diffusion tensor Λ. The gradient scheme is then given by Initialization:

s (0) D (x) = S(x, Π D p (0) (x)), (5.6a)
For n = 0, . . . , N -1 :

u (n+1) -ū D ∈ X D,0 , v (n+1) -v D ∈ X D,0 , (5.6b)
p (n+1) = u (n+1) -v (n+1) , s (n+1)
D (x) = S(x, Π D p (n+1) (x)), (5.6c)
δ (n+ 1 2 ) D s D (x) = s (n+1) D (x) -s (n) D (x) δt (n+ 1 2 ) , (5.6d)
Ω Φ(x)δ (n+ 1 2 ) D s D (x)Π D w(x) +k 1 (x, s (n+1) D (x))Λ(x)∇ D u (n+1) (x) • ∇ D w(x) dx = (5.6e) 1 δt (n+ 1 2 ) t (n+1) t (n) Ω f 1 (x, t)Π D w(x)dxdt, ∀w ∈ X D,0 , Ω -Φ(x)δ (n+ 1 2 ) D s D (x)Π D w(x) +k 2 (x, s (n+1) D (x))Λ(x)∇ D v (n+1) (x) • ∇ D w(x) dx = (5.6f) 1 δt (n+ 1 2 ) t (n+1) t (n) Ω f 2 (x, t)Π D w(x)dxdt, ∀w ∈ X D,0 .
In the case of the Richards equation, the scheme is obtained by replacing one of the discrete conservation equations (5.6e) or (5.6f) by the imposed value for the pressure.

Note that this scheme leads to sparse linear systems if Newton's method is used to solve (5.6) and that this scheme is implementable in practice. Indeed, let us consider the equations corresponding to a test function w given by w i = 1 for some i ∈ I, w i = 0 for all i ∈ I \ {i}, w j = 0 for all j ∈ J. Let us denote by K the support of the basis function χ i . Then the discrete equations thus obtained only depend on the unknowns u

(n+1) i and v (n+1) i
, where i is such that the support of χ i is also K, and on the unknowns u (n+1) j and v

, where j is such that the support of ξ j is some σ, where σ is a face of K. Let us now consider the equations corresponding to a test function w given by w j = 1 for some j ∈ J, w j = 0 for all j ∈ J \ {j}, w i = 0 for all i ∈ I. where j is such that the support of ξ j is a face of K or L, and on the unknowns u

(n+1) i and v (n+1) i
where i is such that the support of χ i is K or L.

Let us remark that the use of the gradient discretization given by (3.1) would not provide an easy discrete problem. Indeed, with this gradient discretization, the discrete gradients involved in (5.6) cannot be reduced to a local equation, using the same computations as that of the proof of Theorem 3.4, since the use of (3.1d) for the test functions gives no way to transform an expression such as Ω k 1 (x, s (x))Λ(x) in (3.1). This would lead to some complications (different gradient discretizations for u and v at least).

The convergence of the scheme (5.6) is proven in [START_REF] Eymard | Gradient schemes for twophase flow in heterogeneous porous media and Richards equation[END_REF] for general gradient discretizations which are coercive, consistent, limit-conforming and compact in the sense recalled in this paper, and which are such that Π D is a piecewise constant reconstruction operator. Therefore, this convergence proof applies to the modified RT 0 scheme (5.6) as stated below. In conclusion, the gradient discretization concept gives a simple way to adapt the mixed finite element context to the above nonlinear problem, by dropping the requirement of the H div conformity for the nonlinear part of the equations.

Theorem 3 . 2 .

 32 Let (T m ) m∈N be a sequence of regular simplicial meshes in the sense of [8, Theorem 3.1 p.14] such that the size h m of the mesh T m , tends to 0 as m → ∞. Let D m = (X Dm,0 , Π Dm , ∇ Dm ) be defined by (3.1) for each m ∈ N, then the family (D m ) m∈N is coercive, consistent, limit-conforming and compact in the sense of the definitions of Section 2.

Theorem 4 . 2 .

 42 to be a norm on X D,0 . The existence and uniqueness of ∇ D u again results from the fact that (4.1c) provides a square linear system, whose solution vanishes if the right-hand-side vanishes. The fact that it defines a norm results, on one hand, from the coercivity property shown in Theorem 4.2 below, and on the other hand, on [18, Proposition 3.1 p. 15], whose consequence is that for given (u i ) i∈I and ∇ D u, there exists one and only one (u j ) j∈J such that (4.1c) hold. Let (T m ) m∈N be a sequence of regular simplicial meshes in the sense of [8, Theorem 3.1 p.14] such that the size h m of the mesh T m , tends to 0 as m → ∞. Let D m = (X Dm,0 , Π Dm , ∇ Dm ) be defined by (4.1) for each m ∈ N, then the family ( D m ) m∈N is coercive, consistent, limit-conforming and compact in the sense of the definitions of Section 2.

Theorem 4 . 3 .

 43 Using the gradient discretization (4.1), the gradient scheme (2.4) for the approximation of Problem (1.1) is equivalent to the Arnold-Brezzi formulation (1.7) of the mixed finite element method.

  Let us denote by σ the support of the basis function ξ j , and by K and L the two cells having the common face σ. Then the discrete equations thus obtained only depend on the unknowns u

  )Λ(x)∇ D u (n+1) (x) • ∇ D w(x)dx in order to get Ω Π D w(x)div(k 1 (x, s (n+1) D (x))Λ(x)∇ D u (n+1) )(x)dx: the H div conformity of k 1 (x, s (n+1) D (x))Λ(x)∇ D u (n+1)is indeed not ensured by the scheme. Therefore the practical implementation of this gradient discretization appears as impossible, if we do not replace Λ by k 1 (x, s (n+1) D

Theorem 5 . 2 (X

 52 Convergence of the numerical scheme). Let ( D m ) m∈N be a sequence of gradient discretizations given by (4.1) with k = 0, and, for allm ∈ N, let (t (n) m ) n=0,...,Nm such that t (0) m = 0 < t (1) . . . < t (Nm) m = T be such that lim m→∞ max n (t (n+1) m -t (n)m ) = 0. Then the scheme (5.6) has at least one solution for any m ∈ N.Moreover, we assume that are given p(0) m , ū Dm , v Dm ∈ X Dm such that the sequences Π Dm p (0) m -p ini L 2 (Ω) , Π Dm ū Dm -ū L 2 (Ω) , Π Dm v Dm -v L 2 (Ω) , ∇ Dm ū Dm -∇ū L 2 (Ω) d , ∇ Dm v Dm -∇v L 2 (Ω) d tend to 0 as m → ∞. Let u m , v m ,s Dm be the piecewise constant time implicit functions equal on (t (n) , t (n+1) ) to u Dm,0 , such that the scheme (5.6) holds for m ∈ N.Then there exist u, v ∈ (L 2 (Ω × (0, T ))) 2 and p = u -v such that, for any t 0 ∈ (0, T ), up to a subsequence,1. (Π Dm u m , Π Dm v m ) converges in L 2 (Ω × (0, t 0 )) 2 to (u, v) as m → ∞, 2. (∇ Dm u m , ∇ Dm v m ) converges in (L 2 (Ω × (0, t 0 )) d ) 2 to (∇u, ∇v) as m → ∞,3. s Dm converges in L 2 (Ω × (0, T )) to s such that s(x, t) = S(x, p(x, t)), for a.e. (x, t) ∈ Ω × (0, T ) as m → ∞, and (u, v) is a weak solution of Problem (5.3).


Hongquan Wang
email: hongquan.wang@insa-rennes.fr

Stéphane Méric
email: stephane.meric@insa-rennes.fr

Sophie Allain
email: sophie.allain@univ-rennes1.fr

Eric Pottier
email: eric.pottier@univ-rennes1.fr

Adaptation of Oh Model for Soil Parameters Retrieval Using Multi-Angular RADARSAT-2 Datasets

Keywords: Multi-angular Measurements, RADARSAT-2 Datasets, Oh Model, Surface Roughness, Soil Moisture

The global context of this study is to improve the robustness of surface geophysical parameter retrieval models. Consequently, the adaptation of the existing retrieval models to radar measurement datasets is very important regarding the improvement of soil parameter retrieval procedures. In this context, the objective of this paper is to empirically adapt the Oh model based on the polarimetric RADARSAT-2 datasets acquired at different incidence angles, along with the simultaneous ground truth. The coefficients in the Oh model are adapted in order to decrease the deviation between Oh model predictions and SAR measurements. Improved agreements are found between the adapted Oh model and backscattering coefficients derived from multiangular (ranged from 24° to 43°) SAR data. This indicates that the model adaptation could provide a way to improve the model behaviours and enhance the geophysical parameter retrieval.

I. INTRODUCTION

Soil surface roughness plays an important role in several physical processes regarding water conservation, soil erosion and surface runoff over agricultural areas. As the conventional soil parameter acquisitions are based on point measurements and considerable time and labour consumption, remote sensing techniques provide an opportunity to overcome these problems with a spatial resolution from centimeter to global scale. Thus, thanks to the signal processing which is called Synthetic Aperture Radar (SAR) processing and applied to radar signals, radar systems can provide images with a very high resolution [START_REF] Cumming | Digital processing of synthetic aperture radar data[END_REF] (RADARSAT-2 has an ultra-high resolution mode of about 3 meters for resolution). Furthermore, compared with optical remote sensing, radar is independent of weather and is able to give subsurface information taking into account the opportunity to use the radar wave polarisation state. Consequently, radar through radar images provides a promising approach to characterize the agricultural soils regarding the sensitivities to soil moisture and surface roughness.

Considering all these assumptions, the soil parameter retrieval can be achieved by considering polarimetric SAR images. Studies on characterization and retrieval of soil parameters like moisture and roughness using polarimetric SAR data have been carried out extensively in the past decades [START_REF] Fung | Microwave Scattering and Emission Models and Their Applications[END_REF]. The SAR signature is affected by not only surface roughness, but also soil moisture. In addition, the radar backscattering signature is dependent on incidence angle, polarization and frequency [START_REF] Ulaby | II -Radar Remote Sensing and Surface Scattering and Emission Theory[END_REF]. Actually, the efficiency of retrieval procedure depends essentially on the scattering model accuracy that describes the soil behavior under electromagnetic wave effects. Thus, the characterization of soil parameters requires the scattering model to be adapted to actual SAR measurements.

One of the most frequently used models is the semi-empirical Oh model [START_REF] Oh | An empirical model and an inversion technique for radar scattering from bare soil surfaces[END_REF][START_REF] Oh | Semi-empirical model of the ensemble-averaged differential mueller matrix for microwave backscattering from bare soil surfaces[END_REF][START_REF] Oh | Quantitative retrieval of soil moisture content and surface roughness from multipolarized radar observations of bare soil surfaces[END_REF], which offers a wide range of applicability criteria. The cross-polarization signature and the ratios of backscattering coefficients are parameterized in terms of surface roughness, soil moisture and incidence angle based on fitting a large number of experimental measurements. However, the discrepancy between the model simulation and acquired backscattering coefficient is frequently large due to the different data scales and limits of model application criteria. Furthermore, when the algorithm is tested over the sites other than those where they are derived from, large deviation also appears due to different soil condition [START_REF] Mozos | Evaluation of empirical and semi-empirical backscattering models for surface moisture estimation[END_REF][START_REF] Merzouki | Evaluation of the Dubois, Oh, and IEM radar backscatter models over agricultural fields using C-band RADARSAT-2 SAR image data[END_REF] that compromises the efficiency of soil parameter retrieval. Thus, adaptation of Oh model is necessary before the implementation of retrieval procedure.

The capability of Oh model to infer surface roughness and soil moisture from RADARSAT-2 datasets is reported in [START_REF] Merzouki | Evaluation of the Dubois, Oh, and IEM radar backscatter models over agricultural fields using C-band RADARSAT-2 SAR image data[END_REF], demonstrating significant disagreements between model predictions and actual data. The potential use of RADARSAT-2 datasets is evaluated in [START_REF] Baghdadi | A potential use for the C-band polarimetric SAR parameters to characterize the soil surface over bare agriculture fields[END_REF][START_REF] Baghdadi | Sensitivity of main polarimetric parameters of multifrequency polarimetric SAR data to soil moisture and surface roughness over bare agricultural soils[END_REF] by examining the sensitivity of polarimetric parameters to surface roughness and soil moisture. The dynamics of linear polarized signature and polarimetric parameters from RADARSAT-2 are evaluated in [START_REF] Adams | Sensitivity of C-band SAR polarimetric variables to unvegetated agricultural fields[END_REF] and the optimal SAR descriptors are determined for agricultural status monitoring. The discrimination between wet and dry moisture conditions is achieved in [START_REF] Bourgeau-Chavez | Assessment of polarimetric SAR data for discrimination between wet versus dry soil moisture conditions[END_REF] using RADARSAT-2 datasets.

The objective of this paper is to evaluate and investigate the potential to adapt the Oh model using polarimetric RADARSAT-2 datasets extended in a multi-angular configuration and acquired in Northwest of France. Section II describes the test site and the multi-angular datasets. The operating characteristics of RADARSAT-2 are introduced and especially the multi-angular datasets are acquired over the study region considering the meteorological conditions. Moreover, the soil moisture is measured by two methods and a conversion relationship is proposed to calibrate the soil moisture measurements. Finally, the surface roughness is measured using a simple chain method and then converted to the classical roughness descriptor. In section III, comparison between model simulation and RADARSAT-2 response is carried out. Then, the adaptation technique is implemented to improve the model performance and to decrease the disagreements of model predictions and multi-angular RADARSAT-2 datasets. Discussions and main conclusions are presented in section IV and section V respectively.

II. DATASET PRESENTATION

A. Multi-angular SAR Data Acquisitions

The studied agricultural areas are located in the site of Pleine-Fougères, near the Mont-St-Michel (N 48°38', E 1°30') in Brittany, France (see Fig. 1). Multi-angular fully polarimetric RADARSAT-2 images (see Fig. 1) were acquired in April 2013 in ascending and descending orbit and in fine quad-polarization mode (Single look complex -SLC products) with spatial resolution of 5.2 m in range and 7.6 m in azimuth and with three different incidence angles of 24°, 31°and 43° (Table 1). The SLC data are firstly extracted as coherence matrix T3 using PolSARpro4.2 [START_REF] Lee | Polarimetric radar imaging from basics to applications[END_REF]. Then, a boxcar filter with a 7-by-7 window [START_REF] Baghdadi | A potential use for the C-band polarimetric SAR parameters to characterize the soil surface over bare agriculture fields[END_REF] is applied to reduce speckle effect. The polarimetric images are ortho-rectified and transformed from slant range geometry to ground range using NEST (ESA SAR toolbox). As the test fields are flat, no compensation of terrain slope is applied to the SAR images. Afterwards, the master SAR image is selected from one of the multi-angular datasets and other images are coregistered to the master image, regardless of descending and ascending orbits. The co-registration average accuracy is achieved around 0.8 pixel. At this step, the processed Pauli color-coded images [START_REF] Lee | Polarimetric radar imaging from basics to applications[END_REF] are analysed in terms of scattering mechanisms. The volume scattering in green is dominant in forested areas and vegetated fields while the double bounce in red occurs predominantly in urban regions. The single bounce in blue is dominant over the expansive agricultural fields. For instance, in Fig. 2, two Pauli images are analysed for different incidence angle datasets and it can be noted that single bounce component tends to increase as the incidence angle increases. In order to locate the bare agricultural fields, the land cover information is collected in collaboration with the COSTEL laboratory from University of Rennes 2 (Climat et Occupation du Sol par TELédétection). The tested bare fields are selected considering their size and their distributions over the common section of multi-angular SAR swaths described in Fig. 1. Collecting the land cover information makes it possible to select 36 bare fields during the April 2013 campaign, the soil moisture measurements by Time Domain Reflectometry (TDR) instrument and the corresponding calibration procedures are presented in the followings. Moreover, the ground truth campaign involving surface roughness measurements by chain and laser methods is also described.

B. Soil Moisture Measurements

The electromagnetic wave traveling speed between two rods installed in TDR instrument depends on the permittivity of the natural soil which is located between these two rods. Hence, by measuring the electromagnetic wave traveling speed in soil mixture, the soil moisture can be derived. Although the TDR acquisition is fast to carry out over the expansive agricultural fields, the TDR measurement accuracy is usually not sufficient to be exploited. Consequently, TDR calibration is often necessary for agricultural soil with organic matters. In our study, the gravimetric method is accomplished as a soil moisture FQ5 FQ11 reference measurement to evaluate the water content of a soil sample [START_REF] Walker | In situ measurement of soil moisture: a comparison of techniques[END_REF]. After measuring the soil moisture status using both TDR and gravimetric method, the calibration relationship can be established based on minimum Root Mean Square Error (RMSE) technique [START_REF] Jaynes | Probability Theory: The Logic of Science[END_REF]. In our case, we investigate the calibration approaches for a common TDR-system: FieldScout-TDR300 (Spectrum Technologies) [START_REF] Spectrum | User's manual: Field Scout TDR 300 soil moisture meter[END_REF]. Therefore, the TDR measured soil moisture is transformed into the calibrated soil moisture. In situ measurements of soil moisture, a depth of 3.8 cm (probe length) over the agricultural fields was measured over the 36 selected bare fields. Furthermore, TDR measurements were implemented over 25 samples distributed homogeneously over each test field so as to represent the general soil moisture status of each field.

C. Surface Roughness Measurements

Surface roughness is commonly described by the standard deviation of the surface height also named Random Roughness Factor s. Several methods for surface roughness measurements have been developed in the past decades [START_REF] Jester | Soil surface roughness measurement -methods, applicability, and surface representation[END_REF] by using i) vertical movable steel needles that make it possible to obtain the relief of the ground surface through a scale grid attached on a board, ii) a laser profiler transmitting pulses and recording the reflected pulses from the soil surface makes it possible to measure the surface roughness with high resolution through the reflected pulse delay. Nevertheless, these two approaches are still really time consuming, which prohibits any use for a mass of agricultural bare fields. Thus, in our study, the chain method proposed by Saleh [START_REF] Saleh | Soil roughness measurements: chain method[END_REF] is adopted as an alternative simpler and faster approach to measure the surface roughness. This method

Pontorson

is based on the principle that as a chain of given length L 1 is placed across a surface, the covered horizontal distance L 2 decreases as surface roughness increases. Therefore, the roughness descriptor called Saleh Roughness Factor (SRF) is defined as: SRF = 100(1-L 2 /L 1 ). As the surface roughness is characterized by s, a relationship is established to transform quantitatively the roughness measurements of SRF to the surface roughness s. Thus, a laser profiler has been developed to assess and calibrate the chain measurements. The effective length of the rail for holding laser instrument is 2.8 m. Several profiles are generated over different surfaces and the local slope effects are removed for each profile. The surface roughness s is computed for each profile as some examples are shown in Fig. 4 for s. Thus, the laser profiler can be used to transform the SRF measurements to the corresponding s values. To achieve this relationship between SRF and s, a regression model [START_REF] Gilley | Random roughness assessment by the pin and chain method[END_REF] is derived as s =a. SRF b where s is given in cm and the coefficients a and b are in function of rainfall amount [START_REF] Gilley | Random roughness assessment by the pin and chain method[END_REF], and also are affected by the linkage length [START_REF] Merrill | Use of chain set for scale-sensitive and erosion-relevant measurements of soil surface roughness[END_REF]. In our study, based on the regression between the laser measurements and synchronous chain sampling data, the coefficients for roughness scale transformation are derived as a= 0.5072 and b = 0.7867. In situ measurements of roughness with a chain of a given length L 1 = 146.5 cm, a 2.2 cm linkage length was measured also over 36 selected bare fields. Over each test parcel, the SRF was measured 10 times uniformly in two perpendicular directions so as to represent the entire field. For the ground campaign in April 2013, all the chain measurements of SRF are described in Fig. 5. Based on the established conversion model, all these measurements can be transformed into random surface roughness descriptor s. The s values range from 1.0 cm to 5.5 cm.

Considering all these in situ measurements, we are consequently able to compare the RADARSAT-2 data with specific retrieval models regarding soil moisture and roughness.

III. ADAPTATION OF OH MODEL TO RADARSAT-2 DATA

A. Comparisons between Oh2004 Model and Radarsat-2 Data

The semi-empirical Oh model [START_REF] Oh | An empirical model and an inversion technique for radar scattering from bare soil surfaces[END_REF][START_REF] Oh | Semi-empirical model of the ensemble-averaged differential mueller matrix for microwave backscattering from bare soil surfaces[END_REF][START_REF] Oh | Quantitative retrieval of soil moisture content and surface roughness from multipolarized radar observations of bare soil surfaces[END_REF] is generally used for soil parameters retrieval. Its last version (Oh2004 model) is based on EM scattering behaviour, using scatterometer measurements and airborne SAR data. The formulation of the crosspolarization signature σ 0 HV and the backscattering coefficient ratios in co-polarized and cross-polarized channels are given by following

With soil moisture mv, surface roughness ks (k = 2 π f /c is the wavenumber), incidence angle θ in radian, g i , m i , n i (i = 1, 2, 3) are the fitting coefficients. The original values of these coefficients, given in Table 2, are derived from the Oh2004 model. In this study, the parameters described through Eqs. ( 1), ( 2) and (3), are forward simulated regarding the satellite configuration characteristics (frequency f, incidence angle θ and polarization) using the mean value (field-level) of measured soil moisture and surface roughness described in section II as inputs. In addition, the co-polarized p and cross-polarized q ratios are computed from the RADARSAT-2 linear polarization channels as p = σ 0 HH / σ 0 VV and q = σ 0 HV / σ 0 VV respectively (σ 0 HV is extracted from the processed coherence matrix as previously described in section II.A). Thus, the results of the simulations are compared with the RADARSAT-2 data derived parameters as shown in Fig. 6, and each point represents the mean value of SAR derived parameters for a given sampling field (it is computed by averaging the values of all the pixels within the target field, except some pixels near the field boundary to remove boundary effects). The significant disagreements between the original Oh2004 model predictions and RADARSAT-2 derived parameters are found in this study. On the one hand, one can see that the linear polarization signature σ 0 HV is overestimated by Oh2004 model (Fig. 6(a)). This overestimation of σ 0 HV is in agreement with the results in [START_REF] Merzouki | Evaluation of the Dubois, Oh, and IEM radar backscatter models over agricultural fields using C-band RADARSAT-2 SAR image data[END_REF]. The average difference between predicted values using original Oh2004 model and measured values from RADARSAT-2 is observed around 3 dB for low incidence angle i.e. 24° and around 2.4 dB for high incidence angle (43°).

On the other hand, the cross-polarized ratio q and the co-polarized ratio p derived from RADARSAT-2 dataset are underestimated using Oh2004 model (Fig. 6(b) and Fig. 6(c) respectively). Firstly, the discrepancy between the original Oh2004 model and extracted RADARSAT-2 q values ranges from 1.6 dB to 3.3 dB in Fig. 6(b). It is also reported in [START_REF] Merzouki | Evaluation of the Dubois, Oh, and IEM radar backscatter models over agricultural fields using C-band RADARSAT-2 SAR image data[END_REF] that the RADARSAT-2 data have larger dispersion from the original Oh2004 model. Secondly, the values of p derived from RADARSAT-2 datasets range from -1 dB to 1 dB, while those simulated by the original Oh2004 model range from -1.3 dB to 0 dB (Fig. 6(c)). The p > 0 dB may due to several reasons, such as the high roughness condition and the weak periodic soil surface structure. This disagreement requires the original model to be adapted. Moreover, the original p model underestimates RADARSAT-2 data, but it is demonstrated in [START_REF] Merzouki | Evaluation of the Dubois, Oh, and IEM radar backscatter models over agricultural fields using C-band RADARSAT-2 SAR image data[END_REF] that p is overestimated by Oh2004 model when ks < 1.5. This different behaviour with [START_REF] Merzouki | Evaluation of the Dubois, Oh, and IEM radar backscatter models over agricultural fields using C-band RADARSAT-2 SAR image data[END_REF] may due to the soil condition difference, as the ground measurements employed in our study cover much wider ranges (for example ground surface roughness) than previous works. The disagreements between RADARSAT-2 datasets and scattering models may be caused by several reasons, such as the bare soil heterogeneity, the vegetation residue influence and the agricultural manipulation disturbance. For example, the spatial variation of soil texture, the high temporal dynamic of soil moisture, the possible soil moisture gradient in term of depth, and the slight topography influence within the sampling fields may partially explain these disagreements. Furthermore, the apparent roughness measured by satellite in the direction orthogonal to flight path is not consistent with the ground truth roughness collection approach. In addition, the multiple scattering within the subsurface also aggravates the model divergence from RADARSAT-2 measurements, dependent on the wave penetration depth. Moreover, our measured surface roughness range may be close to the limit of Oh model validity. Nevertheless, except for these mentioned factors, the error in the ground truth measurement must be taken into account. In this work, the surface roughness is measured initially by the chain method like SRF, and then the measured SRF is transformed to s. As described in section II.A, the datasets of simultaneous chain and laser measurements are not large, which may cause less representativeness of the calibration relationships. Furthermore, a conversion model is used in this study to calibrate the TDR measurements. This conversion model may also cause errors in deriving the corresponding soil moisture. Considering all these points of view, this study proposes to evaluate and to investigate the potentials to adapt the original Oh2004 model in the aim of estimating soil moisture using multi-angular polarimetric RADARSAT-2 datasets.

B. Methodology for the Adaptation

In the Oh2004 model (equations from (1) to ( 3)), the general function (1 -exp (𝑚 ! (𝑘𝑠) ! ! )) is selected as the indicator of surface roughness influence. This function is in agreement with the phenomenon that backscattering signature is towards zero for very flat surface, and becomes independent of surface roughness for very rough condition [START_REF] Oh | Semi-empirical model of the ensemble-averaged differential mueller matrix for microwave backscattering from bare soil surfaces[END_REF]. Consequently, the new idea is to adapt the coefficients in the surface roughness term (g 1 , m 1 , n 1 for σ 0 HV ; g 2 , m 2 , n 2 for q; and g 3 , m 3 , n 3 for p) using iterative Levenberg-Marquardt method [START_REF] Seber | Nonlinear Regression[END_REF] in order to decrease the deviation of Oh2004 model prediction from actual data. The adaptation process has been performed by monitoring the error between the model outputs and the targets of field-level actual values derived from RADARSAT-2 data (considered as field-level targets for the Levenberg-Marquardt). The iterative process is terminated when a minimum of such error was reached. In the adaptation, 50% pixels from the bare agricultural fields are randomly selected and the reminder pixels shall be used in the validation process. The adapted coefficients for crosspolarization signature σ 0 HV , cross-polarized ratio q and co-polarized ratio p are achieved in Table 2. The sign of the coefficients is the same as the original one but the magnitude of some coefficients is changed. The new coefficient values are presented in Table 2 and the final results are shown in Fig. 7. Also, the fitting coefficients g i , m i , n i (i = 1, 2, 3) influences are addressed for each polarimetric measurement.

Cross-polarization σ 0 HV

In the adaptation process of σ 0 HV to RADARSAT-2 data, the angular sensitivity and soil moisture dependence are kept. The adaptation is mainly focused on the surface roughness characterization term (1 -exp (𝑚 ! (𝑘𝑠) ! ! )) and a scale factor g 1 . Therefore, the dependence of σ 0 HV regarding the roughness parameter ks is described in Fig. 8. Moreover, the improvement of the adapted model over the original one is clearly shown in Fig. 8 where the simulation and RADARSAT-2 data have the same scale values. It is also obvious that the disagreement is larger at low incidence angle (Fig. 8(a)) than at high incidence angle (Fig. 8(c)) even after adapting the model. Indeed, these disagreements are induced since ks > 3.5 is beyond the optimal application condition of Oh2004 model (the SAR signature is less sensitive and close to saturation at very high surface roughness condition) whatever the incidence angle value is. Nevertheless, at high incidence angle condition (such as Fig. 8(c)) the cross-polarization signature is more dominated by surface roughness than soil moisture. Furthermore, Fig. 8(a) shows some abnormal parcels, which are below both the original and adapted model. The low cross-polarization backscattering from these parcels under incidence 24° is due to the low soil moisture (mostly around 15% as mentioned in section II.B), implying that σ 0 HV is also predominated by soil moisture at low incidence angle. In the condition of ks < 2.5, several parcels with mv > 15% matched better with original model than with the adapted model (the study in [START_REF] Merzouki | Evaluation of the Dubois, Oh, and IEM radar backscatter models over agricultural fields using C-band RADARSAT-2 SAR image data[END_REF] also found the Oh model performance in terms of RADARSAT-2 data is improved significantly at condition of mv > 25%), but 85% of the parcels are observed more suitable to the adapted model. It is found that the global discrepancy after the model adaptation decreases to 1.5 dB for low incidence angle and to 1.2 dB for high incidence angle. In this sense, the cross-polarization backscattering signature at high incidence angle acquisition shows great potentials to improve the surface roughness characterization. q (dB) from RS2 q(dB) from adapted Oh2004 24

°31 °43 °-2 -1 0 1 2 -2 -1 0 1 2 p (dB) from RS2 p(db) from adapted Oh2004 24 °31 °43 ° (a) (b) (c)
Fig. 8 Original vs adapted Oh2004 model for cross-polarization σ 0 HV using incidence angle (a) 24°; (b) 31°; (c) 43°

Cross-polarized ratio q

The cross-polarized ratio is very sensitive to roughness variation (Fig. 9). In order to decrease the discrepancy between Oh2004 model and RADARSAT-2 data, the m 2 and n 2 coefficients are adapted to be smaller to get the surface roughness description term (1 -exp (𝑚 ! (𝑘𝑠) ! ! )) to become more dynamic to surface roughness variations. The discrepancy of adapted Oh2004 model from RADARSAT-2 data becomes lower when using high incidence angle acquisition rather than using low incidence angle acquisition. This is also in agreement with the observing results that the surface roughness over bare agricultural fields is more suitable to be described using higher incidence angle acquisitions [START_REF] Baghdadi | Potential of ERS and RADARSAT data for surface roughness monitoring over bare agricultural fields: application to catchments in northern France[END_REF]. This incidence angle influence on cross-polarized ratio q is significantly clear in Fig. 9 (and also Fig. 7(b)). It is noted that at low incidence angle 24° (Fig. 9(a)), the adapted model also underestimates the actual data, even the discrepancy becomes less than the original one. At high incidence angle 43° and in the case of ks > 2.8 (Fig. 9(c)), the adapted model overestimates the actual data while the original model underestimates the actual values, but the absolute divergences are similar. However, for both original and adapted Oh2004 model, the general bias of cross-polarized ratio q is larger for smooth surface than for rough surface, in accordance with the Oh model evaluation in [START_REF] Merzouki | Evaluation of the Dubois, Oh, and IEM radar backscatter models over agricultural fields using C-band RADARSAT-2 SAR image data[END_REF]. The systematic bias between adapted q model and RADARSAT-2 data allows us to conclude a dependence on incidence angle.

Co-polarized ratio p

The co-polarized ratio p is less sensitive to roughness variation as shown in Fig. 10. By model adaptation, the m 3 , n 3 values are in accordance with the original Oh2004 model; only the coefficient g 3 is changed to 1.15 (Table 2). This unchanged surface roughness term implies that the dependence of co-polarized ratio p on surface roughness is weak. For small incidence angle 24° (Fig. 10(a)) and moderate incidence angle 31° (Fig. 10(b)), it is noted that when ks < 3.2, the adapted model agrees with the RADARSAT-2 data. Nevertheless, for high incidence angle 43°, both the model before and after adaptation underestimate the RADARSAT-2 datasets, even though the discrepancy decreases for adapted model. When ks > 3.2, the adapted model overestimates the data whatever the incidence angle value is, and the absolute bias is similar to the original model. Thus, even though the general performance of adapted Oh2004 model is better than the original model, it is also with drawbacks since the limitation of SAR derived parameters under high roughness condition (such as ks > 3.2) can not be overcome.

A. Model Validation

The performances of the adaptation results are assessed based on RMSE which is expressed in dB

2 1 1 RMSE= ( ) N i i i O P N = - ∑ (4)
Where the index i describes the field number, O i is the extracted value in dB from RADARSAT-2, P i is the simulated value in dB from models (original and adapted respectively), and N is the total number of fields. The RMSE is able to quantify the magnitude of average modelling error for each model and for each incidence angle. The results are mainly compared in Fig. 11.

It is indicated that RMSE decreases after the Oh model adaptation to RADARSAT-2 datasets. For cross-polarization signature σ 0 HV (Fig. 11(a)), the RMSE at small incidence angle is larger than that at higher incidence angle even after the model adaptation. In the same way, Fig. 11(b) shows the improvement for the cross-polarized coefficient q using multi-angular data, indicating that the RMSE decreases by around 50% owing to the model adaptation. The deviation of q model from actual SAR data at small incidence angle is also larger than at higher angle. Finally, it is shown in Fig. 11(c) that RMSE decreases by around 40% after the p model is adapted. In contrary to the behaviours of σ 0 HV and q, the RMSE of original p model slightly decreases at incidence angle 31° and then increases at 43°. After the p model is adapted, the RMSE increases with incidence angle. The comparison of the model behaviours indicates that the adapted σ 0 HV and q are more suitable to characterize the surface roughness at high incidence angle. This indicates that the sensitivity of polarimetric parameters to surface roughness increases with increasing incidence angle, which is in agreement with the theoretical prediction [START_REF] Baghdadi | Potential of ERS and RADARSAT data for surface roughness monitoring over bare agricultural fields: application to catchments in northern France[END_REF]. Based on examining these RADARSAT-2 datasets, the original cross-polarization σ 0 HV model is found to overestimate RADARSAT-2 data by 2.4 -3 dB. This error is larger than the RMSE error of 2.34 dB in [START_REF] Oh | Semi-empirical model of the ensemble-averaged differential mueller matrix for microwave backscattering from bare soil surfaces[END_REF]. The original q model underestimates the RADARSAT-2 data by around 1.6 -3.3 dB (in [START_REF] Oh | Semi-empirical model of the ensemble-averaged differential mueller matrix for microwave backscattering from bare soil surfaces[END_REF] RMSE equals to 2.06 dB), and original p model underestimates the RADARSAT-2 data by around 0.6 dB (in [START_REF] Oh | Semi-empirical model of the ensemble-averaged differential mueller matrix for microwave backscattering from bare soil surfaces[END_REF] RMSE equals to 0.82 dB). The divergences between original Oh2004 model and RADARSAT-2 data may due to the fact that the original Oh2004 model is derived using datasets from scatterometer and airborne SAR, different from RADARSAT-2 image acquisition system. After adapting the Oh2004 model to RADARSAT-2 datasets, the RMSEs decrease to around • 1.4 dB for cross-polarization σ 0 HV , • 0.8 -1.6 dB for cross-polarized ratio q, • 0.32 -0.45 dB for co-polarized ratio p.

Thus, the adapted Oh model shows less divergence from the RADARSAT-2 data.

As incidence angle increases, the RMSE decreases for both σ 0 HV model and q model (in contrary with p mode). Thus, it is suggested to use σ 0 HV and q formulations at high incidence angle for surface roughness characterization. These findings illuminate the importance of cross-polarization signature for surface roughness quantification over bare agricultural environments. As in C-band, the microwave penetrating ability is weaker than in L-band [START_REF] Nolan | Penetration depth as a DINSAR observable and proxy for soil moisture[END_REF]. Thus, the cross-polarization signature which represents the multiple scattering components may be predominated more by surface roughness than soil moisture.

The adapted Oh2004 model provides a possible application when surface roughness satisfies 1.3 ≤ ks ≤ 5.6 with moderate soil moisture and incidence angle ranging from 24° to 43°. In addition, the adapted Oh2004 model suggests a potential use of multi-angular SAR data for improving the effectivity of backscattering model. Nevertheless, further work is still necessary to test the adapted Oh2004 model using large databases acquired by RADARSAT-2 (and/or other platform) over various study conditions in order to validate the robustness of the adapted Oh2004 model.

V. CONCLUSION

The objective of this study is to evaluate and adapt the original semi-empirical Oh2004 model based on the multi-angular RADARSAT-2 datasets and simultaneous ground truth measurements over bare agricultural soils. The adapted model maintains the sensitivities to the variation of incidence angle and soil moisture. Meanwhile, the sensitivity to surface roughness along with the data scale is determined to minimize the model discrepancy from actual SAR measurements. With the improved surface roughness characterization ability for the bare agricultural fields, the adapted σ 0 HV and q model are proposed for surface roughness retrieval, preferentially using high incidence angle RADARSAT-2 data.

Fig. 1

 1 Fig. 1 RADARSAT-2 multi-angular image swaths over the Pleine-Fougères region

Fig. 2

 2 Fig. 2 Comparisons between two Pauli RGB SAR images and the corresponding optical area for two incidence angles (double bounce (red), volume (green) and single bounce (blue))

  For the ground campaign in April 2013, the soil moisture distributions on 19/04/2013 and 23/04/2014 and the daily variations in April 2013 are shown in Fig. 3. Based on these measurements, the soil moisture ranges can be derived: • From 23% to 32% with standard deviation around 4% on 19/04/2013. • From 15% to 24% with standard deviation around 4% on 23/04/2013. Besides the soil moisture, another important soil parameter is the surface roughness (a) Distribution on 19/04/2013 (b) Distribution on 23/04/2013 (c) Daily variation in April 2013 Fig. 3 The distributions and daily variation of soil moisture measurements in April 2013

  mv (%) 24° incidence angle (23/04) 43° incidence angle (20/04)

  (a) Low surface roughness profile (b) High surface roughness profile Fig. 4 Surface roughness s after slope correction (a) 19/04/2013 (b) 23/04/2013 Fig. 5 The distribution of measured surface roughness SRF in April 2013

Fig. 6

 6 Fig. 6 Comparisons of the RADARSAT-2 datasets with Oh2004 model

Fig. 7

 7 Fig. 7 Comparisons of the RADARSAT-2 datasets with the adapted Oh2004 model

Fig. 9 Fig. 10 Original

 910 Fig.9Original vs adapted Oh2004 model for cross-polarized ratio q using incidence angle (a) 24°; (b) 31°; (c) 43°

  Fig. 11 Comparisons of RMSE (in dB) between original and adapted Oh2004 models

TABLE 1 MULTI

 1

			-ANGULAR RADARSAT-2 DATASETS
	Date (2013) Beam mode Incidence angle (°) Images orbit (Asc/Des)
	23/04	FQ5	24	A
	23/04	FQ11	31	D
	20/04	FQ24	43	D

TABLE 2 ORIGINAL

 2 AND ADAPTED COEFFICIENT OF THE OH MODEL

	Model		σ	0	hv			q			p
	Coefficient	g 1		m 1	n 1	g 2	m 2	n 2	g 3	m 3	n 3
	Original	0.11	-0.32	1.8	0.095	-1.3	0.9	1	-0.4	1.4
	Adapted	0.11	-0.21	1.3	0.17	-0.71	0.75	1.15	-0.4	1.4


Hongquan Wang
email: wang@insa-rennes.fr

Stéphane Méric
email: stephane.meric@insa-rennes.fr

Sophie Allain
email: sophie.allain@univ-rennes1.fr

Eric Pottier
email: eric.pottier@univ-rennes1.fr

Monitoring Bare Agricultural Soil: Comparison between Ground Based SAR PoSAR System Measurements and Multi-angular RADARSAT-2 Datasets

The objective of this study is to evaluate the surface roughness and soil moisture effects on polarimetric parameters under multi-incidence angle condition. In order to achieve this objective, we propose to use Ground Based (GB) SAR experiments, in which the specific levels of surface roughness and soil moisture are set deliberately. The results show that image acquired with high incidence angle is more sensitive to surface roughness than that with low incidence angle. Moreover, the results we obtained with our GB-SAR system (called PoSAR for Pocket SAR) are in agreement with the behavior of multi-angular RADARSAT-2 data that we have previously obtained. This indicates that the combination of various incidence angles may help us to improve surface roughness monitoring accuracy and also to get a better estimation of soil parameters.

Introduction

Polarimetric SAR provides an opportunity to monitor the spatial and temporal distributions of the soil parameters in both regional and continental scale. The SAR response from bare soil is in functions of target parameters (such as roughness and soil moisture) and system parameters (such as frequency, polarization and incidence angle). Theoretical or empirical models are used to convert the measured backscattering coefficients into surface roughness and soil moisture. The multi-angular SAR images have potentials to improve the estimation of soil parameters [START_REF] Zribi | A new empirical model to retrieve soil moisture and roughness from C-band radar data[END_REF]. Multi-angular and multi-polarization configurations were compared in [START_REF] Molineaux | Using multi-angular radar data to discriminate the influence of rough surface scattering on soil moisture inversions over bare soils[END_REF], indicating the sensitivity of multi-angular configuration to surface roughness is ten times than multipolarization, which encourages the continuous exploration of multi-angular inversion models. It was proposed in [START_REF] Zribi | A new empirical model to retrieve soil moisture and roughness from C-band radar data[END_REF] an original surface roughness parameter to include horizontal and vertical roughness statistical information, and found that backscattering difference between two angular SAR acquisitions depends more on this surface roughness indicator than soil moisture. In addition, [START_REF] Rahman | Mapping surface roughness and soil moisture using multi-angle radar imagery without ancillary data[END_REF] suggested to use an additional image acquired under dry condition, which is possible to consider as a unique indicator of surface roughness. As a result, horizontal and vertical roughness can be separated respectively. Nowadays, operational spaceborne and airborne SAR provide major observational images to monitor soil status. Nevertheless, during the time interval of multiangular image acquisitions, the surface status (especially soil moisture) may changes. Moreover, it is quite difficult to accomplish all the ground truth measurements over large areas at the same time. These are obvious constraints with currently orbiting radar systems which can not acquire multi-angular images at the same time. To overcome this drawback, the Ground Based (GB) SAR system provides an opportunity to gap the spaceborne and airborne SAR system limitations. It is capable of acquiring multi-incidence angle images in quasi-invariant soil condition. In this study, we aim to analyze the incidence angle influence on the sensitivity of polarimetric parameters to bare soil. The GB-SAR experiments are analyzed and compared with the behaviors of RADARSAT-2 multi-angular datasets. Section 2 describes the GB-SAR PoSAR experiment and the soil conditions we set up for the measurement. The RADARSAT2 data we used for the comparison are also introduced. In section 3, the angular response of SAR descriptors from the GB-SAR experiment are analyzed and compared with the RADARSAT-2 data obtained under the similar conditions. Finally, main conclusions are presented in section 4.

Experiment descriptions 2.1 GB-SAR system measurement

A very high resolution GB-SAR PoSAR system is developed in the IETR [START_REF] Ferro-Famil | Posar: A vhr tomographic gb-sar system application to snow cover 3-d imaging at x and ku bands[END_REF]. This system is comprised of vector network analyzer, accurate positioning system, antennas, and 3-m rail which undertakes the system. The signal is generated by vector network analyzer in steppedfrequency mode, and then transmitted during system movement which is controlled by the accurate positioning system. The GB-SAR system obtain each line of the image in a stationary state, thus the Doppler centroid is exactly zero. In addition, the value of Doppler rate can be derived by the geometry configuration between the system and target. These accurate estimations of Doppler parameters is essential for signal synthesis. The GB-SAR PoSAR campaign was carried out on 16 October 2013 over the experimental site of Monterfil (N 48 • 4 , E1 • 58 ). The system is mounted on a scaffold with height 6.3 m as shown in Figure 1 We will see in the paragraph 2.2 that these incidence angles correspond to the RADARSAT-2 acquisition incidence angles. Nevertheless, due the limitation of platform height, the incidence angle variation inside each pixel (Figure 2(b)) is much larger than RADARSAT-2 (the incidence angle variation within full swath is small for RS2). This large variation leads an error in angular response analysis. Although this, the GB-SAR PoSAR system operating parameters are set as much as possible to make the measurements be comparable with RADARSAT-2 images.

Surface roughness was measured by a simple and fast chain approach proposed by [START_REF] Saleh | Soil roughness measurements:chain method[END_REF]. It is based on the justification that when a chain of given length L1 (146.5 cm with 2.2 cm linkage length in our case) is placed straightly along a surface, the covering horizontal distance L2 decreases as surface roughness increases. Therefore, the Saleh Roughness Factor (SRF) is defined as: SRF = 100(1 -L2/L1). This SRF roughness descriptor is used in our study. The SRF was measured 10 times uniformly in two perpendicular directions so as to represent the entire field (Figure 3(a)). Moreover, considering the surface roughness is frequently characterized as random roughness factor (s: root mean square height of surface profile) measured using laser/pin roughness meter, a transformation relationship s(cm) = a SRF b is proposed in [START_REF] Gilley | Random roughness assessment by the pin and chain method[END_REF]. The coefficients a and b are in function of rainfall amount [START_REF] Gilley | Random roughness assessment by the pin and chain method[END_REF], and also affected a lot by the linkage length as reported [START_REF] Merrill | Use of chain set for scale-sensitive and erosion-relevant measurements of soil surface roughness[END_REF]. In our study, the coefficients for roughness scale transformation are derived as a = 0.5072, b = 0.7867 based on the regression between our measurements using laser and synchronous sampling data using chain approach (Figure 3 During each GB-SAR PoSAR operation, in situ measurements of soil moisture and surface roughness are implemented, in the similar manner as RADARSAT-2 ground campaign, except for the field dimension here is much smaller.

RADARSAT-2 data acquisition

The study agricultural areas using RADARSAT-2 data are located in the site of Pleine-Fougères, near the Mont-St-Michel (N 48 • 38 , E1 • 30 ) in France. Multi-angular polarimetric RADARSAT-2 images were acquired in April 2013 in fine-quad polarization mode with six incidence angles: 24 • , 31 • , 33 • , 40 • , 43 • and 49 • . The test bare fields are flat and are selected considering their size and their distributions over the common section of multiangular SAR swaths. in situ measurements of surface roughness and soil moisture were carried out in coincidence with the SAR acquisitions over 34 bare fields. As for the GB-SAR campaign, the surface roughness was measured use a simple chain approach and then transformed to conventional roughness descriptor. Moreover, the soil moisture was measured using the same calibrated TDR as in GB-SAR experiments.

SAR data analysis

The results we obtained about the sensitivity of GB-SAR PoSAR data to soil characteristics at different incidence angles are presented. Moreover, the behaviors of GB-SAR PoSAR experiment is compared with the findings in RADARSAT-2 datasets [START_REF] Wang | Soil moisture characterization using multi-angular polarimetric Radarsat-2 datasets[END_REF].

GB-SAR PoSAR data processing

The received signal is synthesized using back projection algorithm [START_REF] Cumming | Digital processing of synthetic aperture radar data[END_REF] and the basic radar scene which is exhibited in the Figure 5(a) is used for calibration (polarimetric and radiometric). As shown in Figure 5(b), the trihedral corner reflector deployed within the scene is well synthesized. Then the polarimetric calibration is implemented on the synthesized images using approach [START_REF] Ainsworth | Orientation angle preserving a posteriori polarimetcalibration SAR calibration[END_REF] to remove the cross talk and channel imbalance. Finally, the radiometric calibration is accomplished by calculating the theoretical radar cross section of trihedral corner reflector [START_REF] Van Zyl | Calibration of polarimetric radar images using only image parameters and trihedral corner reflector responses[END_REF] and comparing with the integral energy of the trihedral within the image.

RADARSAT-2 data processing

RADARSAT-2 data are extracted as coherence matrix T3 using PolSARpro4.2. A boxcar filter with 7 × 7 windows is applied to reduce speckle. The polarimetric images are ortho-rectified using NEST.4B (ESA SAR toolbox). Then, the master image is selected from the multi-angular datasets. Other incidence angle images are considered as the slave images and co-registered to the master image.

The common section of multi-angular image swath is extracted for further analysis and model development. The half of the pixels inside reference parcels are stochastically selected. Thus, the SAR descriptors are statistically calculated for each study parcel.

Sensitivity comparisons between PoSAR and RADARSAT-2

For the multi-angular analysis, the PoSAR image pixels located in the same column are considered to have the same incidence angle. The linear backscattering coefficients in HH polarization are analyzed in Figure 6, implying the rougher surface s=3.4 cm and the original very smooth surface can be separated from other roughness status. In contrary, the moderate surface roughness s=2.5 cm is confused with the smooth roughness s=1.6 cm. It is noted obviously that at high incidence angle condition, the four roughness status are easier to discriminate than at low incidence angle. In the same way, four surface roughness status are selected for the RADARSAT-2 study, considering the roughness levels we used in the GB-SAR PoSAR experiment. The discriminations of four roughness status using RADARSAT-2 data at low incidence angle and high incidence angle are shown in Figure 7, also indicating that the four roughness fields are easier to be separated at high incidence angle. However, we have to notice that the soil moisture is different between the GB-SAR PoSAR experiments and the RADARSAT-2 datasets. Nevertheless, this soil moisture difference does not affect the roughness discrimination ability at high incidence angle (backscattering coefficient is dominated by surface roughness). Thus, the behaviors of GB-SAR PoSAR experiment and RADARSAT-2 data are in accordance with each other, and also agree with the conclusions in [START_REF] Adams | Sensitivity of C-band SAR polarimetric variables to unvegetated agricultural fields[END_REF]. Besides the backscattering coefficient, the polarimetric parameter is also used to discriminate different soil status. Based on the difference between double and volume scattering power, a polarimetric parameter DERD is defined in [START_REF] Allain | Two novel surface model based inversion algorithms using multi-frequency polsar data[END_REF] to characterize the surface roughness. It is shown in Figure 8 the separation among four roughness status using GB-SAR and RADARSAT-2 measurements. The DERD decreases with surface roughness in two cases, and the GB-SAR data is more robust to discriminate the surface roughness.

Conclusions and perspectives

In this study, the sensitivity of GB-SAR measurement to soil characteristics are evaluated under different incidence angles and compared with the findings we obtained using RADARSAT-2 multi-angular datasets. By using the GB-SAR PoSAR system, the multi-angular polarimetric SAR data can be obtained in a quasi-invariant soil condition, which is an indispensable assumption for multi-angular inversion model. Furthermore, different surface roughness and soil moisture condition can be set deliberately, which is superior to field-scale spaceborne SAR conditions. The results we obtained using our GB-SAR PoSAR measurement indicate that at higher incidence angle, the backscattering signature is more dominated by surface roughness. The polarimetric parameter (such as DERD [START_REF] Allain | Two novel surface model based inversion algorithms using multi-frequency polsar data[END_REF]) has the potential to discriminate roughness at high incidence angle. This angular response pattern in GB-SAR PoSAR measurements is in accordance with RADARSAT-2 datasets.

Figure 1 :

 1 Figure 1: PoSAR system deployment on scaffold.

Figure 2 :Figure 3 :

 23 Figure 2: GB-SAR PoSAR experiment (a) incidence angle in the scene; (b) incidence angle variation within a pixel.

  (b)). Based on the original very smooth soil status, three other roughness conditions (s = 1.6, 2.5, 3.4) are set by the tillage operation (Figure 4(a)).

Figure 4 :

 4 Figure 4: Different (a) surface roughness; (b) soil moisture setting up in GB-SAR experiments.

Figure 5 :

 5 Figure 5: (a) The trihedral corner reflectors within scene (ridge length=30cm); (b) HH SAR image after synthesizing signal.

  High θ = 43 •

Figure 6 :

 6 Figure 6: Four surface roughness discriminations based on GB-SAR PoSAR data (soil moisture around 30%)

  High θ = 43 •

Figure 7 :

 7 Figure 7: Corresponding four surface roughness discriminations based on RADARSAT-2 data (soil moisture around 15%).

Figure 8 :

 8 Figure 8: Polarimetric parameter DERD comparison.