Yvon Maday
email: maday@ann.jussieu.fr

Benjamin Stamm
email: stamm@math.berkeley.edu

LOCALLY ADAPTIVE GREEDY APPROXIMATIONS FOR ANISOTROPIC PARAMETER REDUCED BASIS SPACES

Keywords: model reduction, reduced basis method, greedy algorithm

Reduced order models, in particular the reduced basis method, rely on empirically built and problem dependent basis functions that are constructed during an off-line stage. In the online stage, the precomputed problem dependent solution space can then be used in order to reduce the size of the computational problem. For complex problems, the number of basis functions required to guarantee a certain error tolerance can become too large in order to benefit computationally from the model reduction. To overcome this, the present work introduces a framework where local approximation spaces (in parameter space) are used to define the reduced order approximation in order to have explicit control over the on-line cost. This approach also adapts the local approximation spaces to local anisotropic behavior in the parameter space. We present the algorithm and present numerous numerical tests.

1. Introduction. The recent progresses in the numerical simulation of physical phenomena obtained through the combination of robust and accurate approximation methods and faster and larger computational platforms have permitted to investigate more and more complex problems with improved reliability. These progresses have, in turn, lead to new demands for the numerical simulations that are not only used to understand a given state but investigate control and optimization problems. Here it is not only a unique (or very few) instance(s) of a given model that is required but a large number of similar computations of a now parameter-dependent model dynamically varying through the values of the parameter under which the behavior of the solution is sought for.

Faster solution algorithms are often not sufficient to achieve these new demands in many engineering applications and reduced numerical methods have been proposed as surrogates to standard numerical approximations of given mathematical models. Among these methods, the proper orthogonal decomposition (POD) and reduced basis (RB) methods are two main classes of higher-order mathematical techniques that have been widely developed during the last decade. These approaches have gained in generality and reliability [START_REF] Binev | Convergence rates for greedy algorithms in reduced basis methods[END_REF][START_REF] Buffa | A priori convergence of the greedy algorithm for the parametrized reduced basis method[END_REF][START_REF] Machiels | Output bounds for reduced-basis approximations of symmetric positive definite eigenvalue problems[END_REF][START_REF] Prud'homme | Reliable realtime solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF].

The basic idea behind these reduced numerical methods is the notion of small Kolmogorov n-width of the set of all solutions obtained when the parameters varies. Indeed POD or RB methods are based on a two steps strategy : the first step (offline stage) allows to select particular instances of the parameters, for which a very accurate approximation of the solution is computed : the associated solutions constitute the reduced basis, in a second step (the on-line stage) the generic solution (for other instances of the parameter) are approximated by a linear combination of these reduced basis elements. The interest of the approach lies in the fact that, in cases the Kolmogorov n-width is small -and there are many reasons a priori, such as regularity of the solutions as a function of the parameters, and a posteriori, revealed by a large number of simulation -these reduced model approximations require very few degrees of freedom and in some sense are close in spirit to other high order approximation methods like spectral methods, but with an improved efficiency.

It is well known that high order methods generally take advantage of a global approach by using basis functions that have a large support that, combined with important regularity of the solution to be approximated (going some times up to analyticity) allows to get, with very few degrees of freedom, a very good accuracy. For most practical design problems in engineering though, the solution is not analytical and most of the time regularity exists but only locally which precludes the interest of global approximations. This is the reason why, for instance, by breaking the global framework to locally piecewise global approaches, the spectral element methods reveals superiority with respect to plain spectral method: a trade off between locality and globality is generally preferred as is demonstrated in e.g. [START_REF] Cantwell | From h to p efficiently: strategy selection for operator evaluation on hexahedral and tetrahedral elements[END_REF] for approximation in spacial direction by spectral element approximations. The same cause imply the same effects in the parameter directions for RB approximations, this the reason why recently, some works have been devoted to proposing ways of adding locality to the reduced methods yielding parameter subdomain domain refinement.

In this paper we shall focus on the certified RB framework for which, the construction of basis functions during the first stage of the algorithm is, as in most of the current approaches, performed through a greedy strategy based on a posteriori error estimator. This concerns both the Galerkin approximation and the empirical interpolation methods. For either high dimensional parameter spaces or spaces with large measure one may encounter the problem that the size of the reduced basis turns out to be larger than desired. Following the lines drawn above, a first idea in this context has been presented by Eftang, Patera and Rønquist [START_REF] Eftang | An "hp" certified reduced basis method for parametrized elliptic partial differential equations[END_REF] and also Eftang and Stamm [START_REF] Eftang | Parameter multi-domain 'hp' empirical interpolation[END_REF] where the parameter space is decomposed into cells where different reduced basis sets are assembled. This approach presents clear advantages in the size of the matricial system that appears in the on-line solution procedure, and corroborates the natural feeling that, in order to approximate the solution at a given parameter, primarily those solutions in the reduced basis corresponding to parameters that are close to the parameter we are interested in are to be involved in the linear approximation.

A drawback of the current approach [START_REF] Eftang | An "hp" certified reduced basis method for parametrized elliptic partial differential equations[END_REF][START_REF] Eftang | Parameter multi-domain 'hp' empirical interpolation[END_REF] however is that, in two adjacent parameter-subdomains, some of the parameters that are selected may be very close. Due to the difference between the CPU time associated with exploiting a RB method in the on-line stage, and that required for constructing the basis elements in the offline step, this leads to the idea that it might be interesting to be able to use, in one parameter subdomain, the parameters that are used in the adjacent ones. Parameter domain decomposition may not thus be the ultimate approach. Another drawback of all current greedy approach is to be unable to master the size of the discrete system that will be solved in the on-line procedure. These remarks have motivated us to investigate the alternative discussed in this paper.

An outline of the paper is as follows. In Section 2 we present the standard greedy algorithm that is widely used in reduced basis computations and introduce then in Section 3 the concept of local approximation spaces that are able to take into account the local anisotropies in the parameter space. The construction of the local approximation spaces is based on a semi-distance function which is constructed in Section 4 through a learning process of the geometry of the manifold of the solution set. Section 5 introduces the concept of a varying train set that is chosen in a optimal fashion using the possibly anisotropic semi-distance function. Finally, in Section 6 we explain some practical aspects of the on-line implementation and Section 7 presents some numerical examples.

2. The greedy algorithm for reduced basis approximations. Let us first introduce a classical greedy-algorithm to have a common ground to present our ideas. In a general setting, it is assumed that for each parameter value µ in a parameter domain P ⊂ R p , a µ-dependent function v(µ) ∈ W of the variable x in a bounded domain Ω ⊂ R d (e.g. d = 2 or d = 3) can be computed. The space W denotes some functional space of functions defined on Ω. In order to fix the ideas this can be either that v(µ) is a solution to a parameter dependent partial differential equation, an approximation of which can be computed by e.g. classical finite element or spectral method. It can also be that v(µ) is a function f (u(µ)) of some u(µ), the evaluation of u(µ) is easy at some few points in Ω but the general knowledge is required for further purposes. Let

S N = {µ 1 , • • • , µ N }
be a collection of N parameters in P and

W N = span{v(µ 1 ), • • • , v(µ N )}
be the approximation space associated to the set S N . A projection/interpolation operator P N : W → W N is supposed to exists that can be either a parameter dependent Ritz-projection (reduced basis methods in case of a parameter dependent PDE) [START_REF] Prud'homme | Reliable realtime solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF][START_REF] Sen | Natural norm" a posteriori error estimators for reduced basis approximations[END_REF], an interpolation operator (empirical interpolation methods in case v(µ) is a function of some u(µ), e.g. [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Grepl | Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations[END_REF] or simply a L 2 -projection. Further, it is assumed that for each parameter value µ ∈ P an error estimator η(µ; W N ) of the approximation of v(µ) by P N (v(µ)) can be computed (this can be e.g. through a posteriori analysis of the residual, see eg [START_REF] Prud'homme | Reliable realtime solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF]). The case where η(µ; W N ) describes the exact error is not excluded in this general framework.

v(µ) = [u(µ)] 4 or v(µ) = e [u(µ)] )
The following sketch represent a typical greedy algorithm: Let Ξ trial be a chosen finite trial set Ξ trial ⊂ P of representative points.

ClassicalGreedy

1. Choose (possibly randomly) µ 1 ∈ Ξ trial , set S 1 = {µ 1 }, W 1 = span{v(µ 1 )} and N = 1, err = max µ∈Ξtrial η(µ; W 1 ). 2. While err > tol 3.

Find µ N +1 = argmax µ∈Ξtrial η(µ; W N ), err = max µ∈Ξtrial η(µ; W N ).

4.

Compute v(µ N +1 ), set S N +1 = S N ∪ {µ N +1 } and W N +1 = span{W N , v(µ N +1 )}.

5.

Set N := N + 1. 6. End while.

Algorithm: Classical greedy algorithm.

Remark 2.1. The introduction of the finite set Ξ trial is due to practical implementation. Actually the definition of this finite set can evolve during the algorithm. We shall elaborate on this in the frame of our approach latter in the paper in Section 5.

Let Ñ denote the final size of the reduced basis space W Ñ such that the accuracy criterion max µ∈Ξtrial η(µ; W Ñ ) ≤ tol is satisfied. Note that this final size is a consequence of the tolerance that has been chosen, hence, for a prescribed tolerance, we have not the hand on the complexity of the RB approximation method.

3. Local approximations. The new approach pursued in this work is in the same spirit of partitioning the parameter set as in [START_REF] Eftang | An "hp" certified reduced basis method for parametrized elliptic partial differential equations[END_REF][START_REF] Eftang | Parameter multi-domain 'hp' empirical interpolation[END_REF], however (i) we do not impose a clear partition of the parameter space but rather collect a global set of sample points S -preliminary constructed (in an off-line stage) with a given tolerance requirement -(ii) we choose a priori the size of the system we want to solve online by selecting an integer N , then (iii) when a reduced basis approximation is to be computed for a certain given parameter value µ ∈ P, we only use the N basis functions whose parameter values lie in a ball

B µ = { μ ∈ P | d(µ, μ) ≤ r(µ)}, (3.1)
for a given semi-distance function d(•, •). Thus, the radius r(µ) is computed in such a way that there are actually N basis functions in the ball. Therefore the local sample space is defined by

S µ = B µ ∩ S = { μ ∈ S | μ ∈ B µ }
with cardinality equal to N1 . Further the local reduced basis approximation space shall be defined by W µ = span{v( μ) | μ ∈ S µ } and its associated local projection by

P µ : W → W µ .
One basic ingredient of this local reduced basis approximation is thus the metric d(•, •); the radius is adjusted to select the total number of modes N for each value of µ whereas the metric can be chosen isotropic or better accounts for anisotropies in the parameter space. The construction of an adapted local metric is a topic itself and is addressed in the upcoming Section 4.

The basic principle of the Locally Adaptive Greedy Approximation is to first select the number N of basis functions we want to get involved in the on-line procedure (i.e. the size of the matrix system to be solved on-line). A standard greedy algorithm, as described by Section 2, is then performed until N + 1 basis functions are selected. Then, the second stage of the algorithm is performed. The basis functions are no longer searched globally but in the above described local balls. The construction of those requires the data of a metric the construction of which is explained in Section 4 (note that, by default, a uniform metric is possible). This metric is updated at each one of the steps of this second stage. Having such a current metric, only local basis functions are considered in order to enrich the set S N : the locality of the basis functions affects the greedy selection of the basis functions since several basis functions can be chosen per iteration. The following procedure is proposed. At any given iteration, the first sample point, say µ, is chosen such that the estimated error η(µ; W µ ) is maximized over the parameter domain as in the traditional setting except that only a local approximation space W µ is considered to compute the error (-estimator). Next, the maximum estimated error η(µ; W µ ) is searched over all parameter values lying outside of the domain of influence of B µ . By the domain of influence of B µ we mean any µ such that µ ∈ B µ or µ ∈ B µ . Outside of the domain of influence of B µ more basis functions can however be added. This condition assures that at most one sample point is added per ball. This procedure is repeated, and more and more parts of the parameter domain are excluded, until the remaining set is the empty set. Given K global sample points S the following pseudo-code results for one iteration, again for a given subset Ξ trial ⊂ P of representative points :

EnrichmentLoop 1. Set Ξ trial = Ξ trial , err = max µ∈ Ξtrial η(µ; W µ ). 2. While Ξ trial = ∅ and err > tol 3. Set µ K+1 = argmax µ∈ Ξtrial η(µ; W µ ), err = max µ∈ Ξtrial η(µ; W µ ). 4.
Compute v(µ K+1 ) and set S K+1 = S ∪ {µ K+1 }.

5.

Set

Ξ trial := Ξ trial \{µ ∈ Ξ trial | µ ∈ B µ K+1 or µ K+1 ∈ B µ }. 6.
Set K := K + 1. 7. End while.

Algorithm: Basis enrichment loop.

4. Anisotropy and local metric -a general finite difference based approach. In order to reduce at most the number of elements involved in the reduced basis approximation, one has to estimate, on the fly, the way the (unknown) solution v(µ) depends on µ.

Indeed, for instance assuming that the parameter domain P is a subset of R 2 , so that µ = (µ 1 , µ 2 ), and assuming that there exists a function ϕ of one variable, such that v(µ) = ϕ(µ 1 + µ 2 ), then the optimal selection of parameters should be sought e.g. along the line µ 1 = µ 2 involving only a one dimensional parameter set. Of course the fact that v depends only on µ 1 + µ 2 is generally not known. What we have to figure out from prior computations, during the greedy algorithm, is if in the solution dependency, there is such an anisotropy. This is a quite classical quest in numerical analysis, and the knowledge of this valuable information implies a large computational reduction. As in other instances, e.g. adaptive approaches like the one used in finite element approximations (see e.g. [START_REF] Danaila | A finite element method with mesh adaptivity for computing vortex states in fast-rotating bose-einstein condensates[END_REF]), the goal is the construction of a (semi-) metric d(•, •) that accounts for the anisotropic behavior of changes of v(µ) with respect to variations in µ. Similarly as in other approaches, the metric is derived from the knowledge of the Hessian of the function v with respect to the parameter µ.

We propose here a general framework to obtain an approximation of the Hessian by finite differences. We always assume that Ξ trial is such that its convex hull spans the parameter space P.

1 Definition of the Hessian. The goal is to define a Hessian matrix H(µ) for each point µ ∈ Ξ trial upon which the metric will be based on. Since the Hessian is based on the reduced basis approximation, it is updated/constructed at each iteration in the algorithm. To do so, compute a reduced basis approximation v rb (µ) of v(µ) based on the current reduced basis, on the stencil µ(α) = µ + p i=1 α i δµ i , with µ ∈ Ξ trial and where δµ i , i = 1, .., p, is a positive small increment in the i th direction in the parameter and α i = -1, 0, 1, such that, at most, two of them are non zero. From these 3 p approximations v rb (µ(α)) of v(µ(α)) on the surrounding stencil µ(α) of µ, one can, by finite differences, evaluate an approximation of the Hessian of v(µ) as a function of µ by taking e.g.

H ii (•; µ) = D ij v rb (µ) on Ω where the operator D ij is defined by D ij v rb (µ) = v rb (µ + δµ i + δµ j ) -v rb (µ -δµ i + δµ j ) -v rb (µ + δµ i -δµ j ) + v rb (µ -δµ i -δµ j ) 4δµ i δµ j .
Remembering that v rb (µ) is a function on Ω, the above formula thus defines a Hessian matrix for each point x ∈ Ω and parameter value µ ∈ P. In order to obtain one Hessian H ij (µ) per parameter point µ ∈ P one needs to define an average operator over the space Ω. If {ϕ n } N n=1 denotes the basis of W µ used to build v rb such that

v rb (x; µ) = N n=1 α n (µ)ϕ n (x),
we propose to define an averaged Hessian over the Ω-space as follows

H ij (µ) = N n=1 α n (µ) D ij α n (µ).
The term D ij α n (µ) describes the anisotropy of the n-th mode while α n (µ) indicates the its weight in the approximation. 2 On-line construction of the metric. Defining a mesh consisting of P -dimensional simplices with vertices being the set Ξ trial allows to interpolate/represent quantities only defined on Ξ trial by means of P 1-finite element functions on the mesh for any parameter value µ ∈ P, since the convex hull of Ξ trial coincides with the parameter space P. We can therefore easily access the interpolated Hessian H(µ), for any µ ∈ P, by means of its interpolated component functions

H ij (µ).
Then, the construction of d(µ 1 , µ 2 ) is based on the Hessian at the points µ 1 and µ 2 . Indeed, we will consider an average of the two local Hessian H(µ i ), i = 1, 2, to construct

H µ1µ2 = H(µ 1 ) H(µ 1 ) F + H(µ 2 ) H(µ 2 ) F ,
where • F denotes the Frobenius norm. The choice of this linear combination between H(µ 1 ) and H(µ 2 ) is motivated by the following argument. The Hessian is derived from a finite difference stencil and the order of the coefficients is not known and might vary in parameter space. In order to take equal quantities, we propose therefore to normalize the Hessian. This does not result in a loss of information since all we want to extract from the Hessian is the local anisotropy, but averaged properly.

One can then perform an eigenvalue decomposition

H µ1µ2 = V Λ V T
where V is an orthogonal and (Λ) ii = λ i a diagonal matrix consisting of the eigenvalues λ i . Consider the diagonal matrix

|Λ| ii = |λ i | λ 2 1 + . . . + λ 2 P , i = 1, . . . , P,
and the associated symmetric positive definite matrix M µ1µ2 = V |Λ| V T to define the semi-distance between µ 1 and µ 2 by

d(µ 1 , µ 2 ) = (µ 1 -µ 2 ) T M µ1µ2 (µ 1 -µ 2 ).
5. Varying train sets based on Hessian. The goal of this section is to present an approach that further allows to reduce the offline computational cost. Here, we focus on the trial set Ξ trial . It is aimed to keep its cardinality as small as possible, but large enough to capture the local geometry of the parametrized system. Since not only the cardinality of Ξ trial but also the locations of its points matters we propose to use the above constructed metric to design problem adapted training sets of appropriate size which will increase with increasing number of basis functions selected.

Firstly, we construct at each iteration of the presented algorithm a new trial set Ξ trial . The cardinality of the Ξ trial being variable and an increasing function of the inverse of the actual error (-estimation) err. Let Q m be the minimal cardinality of Ξ trial at the beginning of the algorithm and Q M the maximal target cardinality at the end of the algorithm once the tolerance criteria err < tol is reached, then we define

Q(err) = Q M -Q m log(tol) log(err) + Q m
to be the number of points at the next iteration. Having the number of elements in Ξ trial for the next iteration specified it remains to construct an updated pointset Ξ trial . The trial set Ξ trial is constructed such that edges of the corresponding unique Delaunay triangulation are uniform in the slightly modified metric

d(µ 1 , µ 2 ) = d(µ 1 , µ 2 ) r(µ 12 ) ,
where

µ 12 = 1 2 (µ 1 + µ 2 ). Taking the metric d(•, •) instead of d(•,
•) is motivated by the fact that a ball defined by (3.1) corresponds (approximately) to a unit ball in d(•, •). This modification of the metric allows to compare unit balls at a global level.

Remark 5.1. Enriching the trial set during the sampling process is not a new idea. An adaptively enriching greedy sampling strategy was proposed in [START_REF] Haasdonk | A training set and multiple bases generation approach for parameterized model reduction based on adaptive grids in parameter space[END_REF]. We however want to mention the featurse that distinguishes our approach from the previously proposed version on a conceptual level. First our approach not only detects where the manifold of the solutions (as a function of the parameters) is complex but also the directions that are important to sample and those that are not. Second in our approach the size of the discrete system to be solved on-line is determined by the user leading to a controlled simulation cost.

6. On-line implementation. It is well known that the stable implementation of the reduced basis method requires the orthonormalization of the snapshots v(µ), µ ∈ S [START_REF] Prud'homme | Reliable realtime solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF][START_REF] Veroy | Certified real-time solution of the parametrized steady incompressible Navier-Stokes equations: Rigorous reduced-basis a posteriori error bounds[END_REF]. In the standard approach, this orthonormalization is prepared off-line and the vectorial space W spanned by the series {v( μ) | μ ∈ S} is not affected by this change of basis. In our alternative method, there are a large number of approximation spaces, almost one for each µ : W µ = span{v( μ) | μ ∈ S µ }. We cannot compute all various orthonormalized basis sets for all possible W µ off-line. This would be too costly and using way too much memory.

Nevertheless, we can still prepare the orthonormalization process so that, the online orthonormalization is feasible. For any µ and µ ′ in S, let us compute off-line the scalar products

ℓ µ,µ ′ = v(µ), v(µ ′ ) (6.1)
where •, • represents an appropriate scalar product : say either L 2 or H 1 type. Note that we can even compute only those scalar products corresponding to pairs of parameters (µ, µ ′ ) that are close, indeed if µ and µ ′ are distant in P, they will never be in the same S μ.

In the on-line stage, once the approximation in W µ will be required, we can easily orthonormalize the basis set {v( μ) | μ ∈ S µ } first by declaring an order in the parameters

{ μ ∈ S µ } = {µ 1 , µ 2 , ..., µ N } (6.2)
then perform a classical Gramm-Schmidt orthonormalization process

ζ 1 = β 1,1 v(µ 1 ), ζ 2 = β 2,2 v(µ 2 ) + β 2,1 ζ 1 , . . . ζ N = β N,N v(µ N ) + N -1 i=1 β N,i ζ i , (6.3)
where the coefficients β n,i are chosen so that ζ n is orthogonal to

span{v(µ i ) | i = 1, .., n -1} = span{ζ i | i = 1, .., n -1}
and the coefficients β n,n are chosen so that ζ n is norm 1. It is well known that these coefficients exist, are unique and that, in order to compute them, we have to know each scalar product v(µ n ), ζ i , i = 1, .., n -1. It should be noticed that computing those coefficients v(µ n ), ζ i on-line depends on the length of the basis vectors representing the basis functions v(µ n ), which should be avoided during the on-line process. Instead, these scalar products can be easily computed from the precomputed coefficients ℓ µi,µj . First, notice that the basis functions {ζ n } N n=1 can alternatively be expressed by the following linear transformation

ζ 1 = γ 1,1 v(µ 1 ), ζ 2 = γ 2,2 v(µ 2 ) + γ 2,1 v(µ 1 ), . . . ζ N = γ N,N v(µ N ) + N -1 i=1 γ N,i v(µ i ), (6.4)
from the set of basis functions {v(µ n )} N n=1 . It is our intention to construct the coefficients γ n,i representing the change of basis functions from {v(µ n )} N n=1 to {ζ n } N n=1 as defined by (6.4) in a stable way and without any operation depending on the length of the vectors representing v(µ n ) to guarantee stability and good on-line performance.

Remark 6.1. Note that the coefficients γ n,i could be obtained based on (6.4) directly. Indeed, each ζ n is orthogonal to span{v(µ i ) | i = 1, .., n -1} allowing to compute the coefficients γ n,i , i = 1, . . . , n -1. This however involves solving a linear systems based on the mass matrix with coefficients ℓ µi,µj = v(µ i ), v(µ j ) , i, j = 1, .., n -1, which can be heavily linearly dependent leading to an unstable scheme.

The coefficients β n,k can be obtained by the following formula mimicking the Gram-Schmidt procedure (6.3) but using the precomputed coefficients ℓ n,k . For the construction of the coefficients β n,j , we proceed as follows (it helps to have in mind that lj,i stands for ζ j , v(µ i ) ) : for all n = 1, . . . , N apply recursively First set

βn,n = 1, βn,k = -lk,n , k = 1, . . . , n -1, so that v(µ n ) - n-1 k=1 βn,k ζ k is orthogonal to span{ζ i ) | i = 1, .., n -1}. Then α = ℓ n,n - n-1 k=1 ( lk,n ) 2 1 2 , β n,k = βn,k /α, k = 1, . . . , n -1, ln,k = β n,n ℓ n,k + n-1 j=1 β n,j lj,k , k = 1, . . . , N.
Note that the construction of α is not subject to instabilities since the sum involves only positive values. In theory, the set of coefficients {γ n,i } N n,i=1 can be obtained by O(N 3 ) operations from the coefficients β n,j , through a triangular process. Indeed, the coefficients γ n,i can be obtained by the following recursive formula

γ n,i =    n-1 k=1 β n,k γ k,i if i < n, β n,n if i = n, 0 if i > n,
for all n = 1, . . . , N.

In practice however, the set of coefficients {γ n,i } N n,i=1 is even not required, as explained in the following.

As mentioned above the set of basis functions {ζ n } N n=1 leads to well-conditioned matrices, the precomputed quantities are nevertheless expressed in the basis {v(µ n )} N n=1 . As an illustration assume that the reduced basis problem to be solved is about the approximation of the solution to a parametrized PDE, written in a variational form as : for a given parameter value µ, find u such that

a(u, v; µ) = f (v; µ)
where it is assumed, for the sake of simplicity that a and f are affine decomposable

a(w, v; µ) = P p=1 g p (µ)a p (w, v), f (v; µ) Q q=1 h q (µ)f q (v).
The ultimate goal is to construct the well-conditioned stiffness matrix a(ζ i , ζ j ; µ) derived from the off-line pre-computation of the series α p i,j = a p (v(µ i ), v(µ j )), for 1 ≤ i, j ≤ N . To do so, we first compute the series αp

i,j = a p (v(µ i ), ζ j ), for 1 ≤ i ≤ N , recursively over j, 1 ≤ j ≤ N : αp i,j = a p (v(µ i ), ζ j ) = β j,j α p i,j + j-1 k=1 β j,k αp i,k .
Then, compute the series a p (ζ i , ζ j ), for 1 ≤ i ≤ N , recursively over j, 1 ≤ j ≤ N : In a first step, Section 7.1, we only test the feature of local anisotropic approximation balls without adapting the train set using the metric. In a second step, Section 7.2, we test the entire algorithm.

a p (ζ i , ζ j ) = β j,j αp i,j + j-1 k=1 β j,k a p (ζ k , ζ j ), (6.5
As a reduced basis model W, we use a parametrized family of functions combined with the L 2 -projection P . The error is computed in an exact manner and no a posteriori estimation is used. For all numerical tests, Ω and P are discretized by a regular lattice of 75 × 75 points and the target tolerance is set to 10 -3 if not otherwise mentioned.

7.1. Numerical results with fixed train set. 7.1.1. Test 1. We start with presenting a numerical example to illustrate the benefit of the local anisotropic approximation spaces. Consider the function

f 1 (x; µ) = exp -(x 1 -0.1(µ 1 -µ 2 )) 2 0.01 - (x 2 -(µ 1 + µ 2 )) 2 0.01 , x ∈ Ω = (-1, 1) 2 , µ ∈ P = [-0.5, 0.5] 2
that exhibits a constant anisotropy of parameters over the whole parameter space. The dependency in the (µ 1 + µ 2 )-direction is ten times stronger than in the (µ 1µ 2 )direction. Figure 7.1 (left) illustrates the evolution of the number of basis functions that are necessary to be computed versus the achieved accuracy in the L ∞ -norm for N = 20 for the anisotropic approach and the results of the same algorithm but using isotropic local approximation spaces as comparison. One can observe that using the anisotropic approach is beneficial in terms of the number of truth solutions to be computed at the off-line stage. Figure 7.1 (right) presents the same quantity but for varying N = 20, 30, 40. Not surprisingly the number of needed basis functions decreases while increasing N and more and more the exponential convergence establishes as for the classical Greedy-algorithm is expected. We can also observe that for a lower N more basis functions can be included per iteration. Figure 7.2 illustrates the local approximation spaces, the radius as a function of the parameter and the sample set S for the anisotropic and the isotropic approach.

Test 2.

The previous example is in some sense idealized since the anisotropy is not varying throughout the parameter domain. In this regard, the next example

f 2 (x; µ) = exp - (x 1 -(µ 2 1 + µ 2 2 )) 2 0.01 - (x 2 -(µ 2 1 + µ 2 2 )) 2 0.01 , x ∈ Ω = (-1, 1) 2 , µ ∈ P = [-0.5, 0.5] 2
is more interesting. It presents a family of parametrized functions where the functions (as functions of x) are constant along concentric circles around the origin in parameter space. As an example, all four corners of the parameter space define the same function.

We compare again the number of required truth solutions for the anisotropic approach compared with the presented algorithm but using isotropic approximation spaces, which is presented in Figure 7.3 (left) for N = 5. We observe that 60% of the number of truth solutions to be computed can be saved. The corresponding local approximation spaces, radii and sample points S are illustrated in Figure 7.4. Figure 7.3 (right) also presents the number of required truth solutions for N = 5 but for different tolerance levels. We observe that starting from a critical value N c of computed truth solutions, that depends on the tolerance, the curves show different behavior. During the iterations before N c the complete parameter region has to be scanned because the tolerance is achieved nowhere in the parameter space. During the iterations after N c however the tolerance is achieved in an increasing region of the parameter space and this part needs no longer to be not scanned. The enrichment of new truth solutions is then limited in space and the final phase of the greedy algorithm converges faster.

Additionally we present in Figure 7.5 the local approximation spaces for the different tolerance levels. One can clearly see how the approximation spaces take into account the local geometry that also allow non-convex approximation spaces.

f 3 (x; µ) = exp - (x 1 -(µ 1 + 3µ 2 )) 2 0.1 + 5|µ 1 + 3µ 2 | - (x 2 -(3µ 1 -µ 2 )) 2 0.1 + 5|3µ 1 -µ 2 | , x ∈ Ω = (-1, 1) 2 , µ ∈ P = [-0.5, 0.5] 2
is interesting in the sense that it presents an almost singularity in parameter space at the origin. The performance of the anisotropic approach, compared to the isotropic version, is presented in Figure 7.6. Again, the anisotropic approaches outperform the isotropic one and about a 56% of computations of truth solutions can be saved. The local approximation spaces, the radius and the sample points with N = 10 is presented in Figure 7.7. We observe a posteriori that the trial space Ξ trial was not sampled fine enough. Indeed in the region around the origin (and the cross for the isotropic version), every trial point is included in the set of sample points. This also explains the sudden drop of the convergence in Figure 7.6, in particular for the isotropic version. This is a known problematic in greedy methods and our approach with adapted/moving trial sets presents a solution to this problem. The corresponding numerical results are presented in the next section. Further we recognize that the scheme detects the cross where the behavior of the function is most singular as can be seen by the chosen sample points in Figure 7.8 (right). Additionally, we present the sample points that were chosen for different functions of the type

f 3,ξ1,ξ2 (x; µ) = exp - (x 1 -ξ 1 (µ)) 2 0.1 + 5|ξ 1 | - (x 2 -ξ 2 (µ)) 2 0.1 + 5|ξ 2 | x ∈ Ω = (-1, 1) 2 , µ ∈ P = [-0.5, 0.5] 2
where ξ 1 and ξ 2 are functions of µ = (µ 1 , µ 2 ). We can observe the adaptive nature of the algorithm. Again, we plot the local approximation spaces for different values of the tolerance levels in Figure 7.9. We observe that the originally non-connected approximation spaces are becoming connected with increased tolerance. 7.2. Numerical results with adapted train sets. In order to compare the accuracy of this version we introduce a fixed test set of sample points Ξ test . In the following examples it consists of a lattice of 75 × 75 points. It corresponds to the fixed trial set of the tests of the previous section. We use FreeFem [START_REF] Hecht | FreeFEM++ Manual[END_REF] in order to adapt the sample points by constructing the desired uniform mesh in the new metric.

It is not expected that the new approach with an adapting trial set uses fewer computed solutions than the approach using a fixed trial set since a non-adapted trial set may leave large errors in parameters that are outside the fixed trial set. On the contrary the adapted set tracks these forgotten parameter values. The benefit will be on the computational efficiency using the framework of a small trial set in the beginning that is gradually increasing. The points of the trials set are however chosen wisely using the information available from the local geometry of the system. Figure 7.10 illustrates the total number of basis functions that are necessary to be computed versus the achieved accuracy of the proposed algorithm with adapting trial sets for the test cases 1, 2 and 3. In each plot, we present the error on the moving trial set, the error on the fixed testing set Ξ test and the accuracy of the version with a fixed trial set (which equals Ξ test ) as comparison.

7.2.1. Test 1. From Figure 7.10 (left), one can observe that the error on the trial set and on the test set of the approach with an adapted trial set are similar. This illustrates that the accuracy is not only guaranteed on the adaptive trial set but also satisfied on the test grid. Both errors are also similar to the error of the algorithm using a fixed trial set as presented in Section 7.1.

The adaptive version is as good as the version on a fixed trial set, but uses less evaluation of the error estimate as is shown in the following table: 7.2.2. Test 2. Next, we consider the second example of the previous section. Figure 7.10 (middle) compares this version (with adapted trial sets) with the version of the previous section using a fixed trial set. One can observe again that the performance is similar than the version using a fixed trial set, and that the accuracy is also satisfied on the fixed test set Ξ test .

Comparing again with the version using a fixed trial set, the version with an adaptive trial set requires less evaluation of the error estimate as is shown in the following table :   Number of trial point evaluations during greedy using adapted trial sets: 37'444 Number of trial point evaluations during greedy using fixed trial sets:

123'750

We observe a gain of 70% fewer error evaluations. Figure 7.12 presents the trial set, local approximation spaces, radii and sample points of this calculation. 7.2.3. Test 3. Finally, we apply the new algorithm with adapted trial set to the third example. As already mentioned, the solution of this example has a singularity and the version with a fixed train set did not resolve the singularity since the train set was not sampled fine enough. We expect the version with the problem-adapted train set to completely resolve the singularity, and thus in contrast using more sample points than the version with a fixed trial set. Figure 7.10 (right) presents the evolution of the number of sample points versus the accuracy of the algorithm. One can observe that this approach (with moving trial set) uses indeed more sample points than the version on a fixed train set as explained above. Finally, Figure 7.13 presents the trial set, local approximation spaces, radii and sample points. We observe that the trial points are in accordance with the sample points, and that a more dense sampling is indeed required around the origin. We observe that the local approximation spaces are now connected also for the tolerance of 10 -3 .

In order to illustrate the benefit of using adaptive trial sets also in this case we consider a test sample of 75 × 75 uniformly distributed points in the region [-0.05, 0.05] 2 around the origin. Figure 7.14 illustrates the error distribution using the online procedure generated using a fixed (left) and an adaptive (right) trial set. The maximum error is 0.043 resp. 0.00146. While the error tolerance is almost satisfied in the latter case, it is clearly not the case for the former approach. Further, the number of error evaluations is given in the following table: Number of trial point evaluations during greedy using adapted trial sets: 215'888 Number of trial point evaluations during greedy using fixed trial sets: 241'875 Thus, the adaptive version uses still less error estimator evaluations and is more accurate in the region around the origin. Also, the error is more equally distributed. 8. Conclusions. We presented a generalization of the classical greedy algorithm that is used in the framework of reduced basis methods and the empirical interpolation method. The presented algorithm introduces local approximation spaces (in the parameter space) that also account for local anisotropic behavior instead of a global approach. The key idea was to consider the N closest basis functions (for a fixed N ) where the semi-distance is measured in an empirically built metric which can be constructed on the fly and which is problem-dependent. Numerical tests illustrate the benefit of this approach.

  Comparison under variation of N .

Fig. 7 . 1 :

 71 Fig. 7.1: Test 1: Accuracy with respect to the number of truth solutions to be computed in comparison with the isotropic approach for N = 20 (left) and the for different values of N = 20, 30, 40 (right).

Fig. 7 . 2 :

 72 Fig. 7.2: Test 1: Local approximation spaces for selected parameter values (left), radius as a function of the parameters (middle) and sample points (right) for anisotropic approach (top) and the isotropic version (bottom, as comparison) with N = 20.

Fig. 7 . 3 :

 73 Fig. 7.3: Test 2: Accuracy with respect to the number of truth solutions to be computed in comparison with the isotropic approach (left) and the for different end tolerances (right). In both cases N = 5.

Fig. 7 . 4 : 2 Tol=10 - 5 Fig. 7 . 5 :

 742575 Fig. 7.4: Test 2: Local approximation spaces for selected parameter values (left), radius as a function of the parameters (middle) and sample points (right) for the presented approach (top) in comparison wit the isotropic version (bottom) for N = 5.

Fig. 7 . 6 :

 76 Fig. 7.6: Test 3: Accuracy with respect to the number of truth solutions to be computed in comparison with the isotropic approach for N = 10.

Fig. 7 . 7 : 2 ξ 1 = µ 1 +3µ 2 , ξ 2 = 1 2 µ 1 - µ 2 - 2 ξ 1 = µ 1 + µ 2 , ξ 2 = µ 1 -µ 2 Fig. 7 . 8 :

 77212221122278 Fig. 7.7: Test 3: Local approximation spaces for selected parameter values (left), radius as a function of the parameters (middle) and sample points (right) for the presented approach (top) in comparison wit the isotropic version (bottom) for N = 10.

Fig. 7 . 12 :

 712 Fig. 7.12: Test 2: Trial points, local approximation spaces for some specific parameter values, radii and sample points for greedy algorithm with adapted trial set and N = 5.

Fig. 7 . 13 :

 713 Fig. 7.13: Test 3: Trial points, local approximation spaces for some specific parameter values, radii and sample points for greedy algorithm with adapted trial set and N = 10.

Fig. 7 . 14 :

 714 Fig. 7.14: Test 3: Error distribution on a region around the origin of the greedy version using a fixed (left) and an adaptive (right) trial set.

  N 2 ) operations. Similarly, the evaluation of the each component of the vector f (ζ i , µ), 1 ≤ i ≤ N can be performed in O(P N ) operations, and the inversion of the associated matricial problem is done in O(N 2 ) further operations.7. Numerical results. Before we start with the numerical tests, an abstract description of the entire proposed local greedy algorithm is given in the following box:

	before finally assembling the entries of the stiffness matrix
	P
	a(ζ i , ζ j ; µ) =	g p (µ)a p (ζ i , ζ j ),	(6.6)
	p=1
	that involves O(P LocallyAdaptiveGreedy
			)

========= Stage 1 ========= 1. Perform a classical greedy algorithm to select N + 1 basis functions. ========= Stage 2 ========= 2. Compute the error estimate η(µ, W µ ) at each point µ ∈ Ξ trial . 3. Compute the metric function d. 4. Enrich the set of basis functions. 5. Create a new trial set Ξ trial . 6. Go to 2. until tolerance tol is achieved. Algorithm: Locally adaptive greedy algorithm.

Note that N has not the same meaning as in the previous section any more; hence the cardinal of S is not N but a K > N generally much larger.

Acknowledgement. The authors want to thank the help and discussions with Frederic Hecht about the adaptation of FreeFem++ to their purpose. This work was also initiated out of the discussions during the WORKSHOP JLL-SMP: "REDUCED BASIS METHODS IN HIGH DIMENSIONS" held in June 2011, organized with the support of the Fondation de Sciences Mathématiques de Paris.

Number of trial point evaluations for the greedy using adapted trial sets: 133'866 Number of trial point evaluations for the greedy using fixed trial sets: 382'500

We observe a gain of 65% fewer error evaluations. Figure 7.11 presents the trial set, local approximation spaces, radii and sample points for this calculation.


Gilles Boulet
email: gilles.boulet@cesbio.cnes.fr

Albert Olioso

Eric Ceschia

Olivier Marloie

Benoît Coudert

Vincent Rivalland

Jonas Chirouze

Ghani Chehbouni

An empirical expression to relate aerodynamic and surface temperatures for use within single-source energy balance models

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

INTRODUCTION

Assessing the turbulent fluxes of latent and sensible heat at the land surface is a crucial issue for both water resource management (computation of evapotranspiration) and meteorological forecasting (evolution of the planetary boundary layer). In order to compute these fluxes at a suitable spatial scale, estimation methods based on the use of remote sensing data are favoured. There is a large array of evapotranspiration estimation methods that use as input or constraint remotely sensed variables such as the NDVI [START_REF] Courault | Review on estimation of evapotranspiration from remote sensing data: From empirical to numerical modeling approaches[END_REF][START_REF] Gowda | ET mapping for agricultural water management: present status and challenges[END_REF][START_REF] Olioso | Future directions for advanced evapotranspiration modeling: Assimilation of remote sensing data into crop simulation models and SVAT models[END_REF]. Evapotranspiration in potential conditions (i.e. not water limited) can be assessed with relatively good precision from in-situ [START_REF] Cleugh | Regional evaporation estimates from flux tower and MODIS satellite data[END_REF] or remote sensing derived [START_REF] Venturini | Estimation of evaporative fraction and evapotranspiration from MODIS products using a complementary based model[END_REF] meteorological data, and, very often, NDVI. But when water stress occurs the latent and sensible heat fluxes are more difficult to assess. In those cases, there is a tight coupling between the evapotranspiration and the radiative surface temperature and, consequently, methods based on remote sensing data in the Thermal Infra Red (TIR) domain are favoured [START_REF] Kalma | Estimating Land Surface Evaporation: A Review of Methods Using Remotely Sensed Surface Temperature Data[END_REF]. Those methods often compute the instantaneous latent heat flux as the residual of the energy budget and, in most cases, an expression of each individual term of the energy budget is proposed [START_REF] Boulet | Monitoring water stress using time series of observed to unstressed surface temperature difference[END_REF]:

LE H G Rn    (1) n.3
While net radiation Rn and, to a lesser extent, soil heat flux G can be expressed directly as a function of the radiative surface temperature, the turbulent fluxes H (sensible heat) and LE (latent heat) depend on a mixed-surface (soil, air and vegetation) temperature source T aero : (2)

where  is the air density, c p the specific heat of air at constant pressure, T air is the air temperature at a reference level above the canopy and r a,h the aerodynamic resistance for heat exchange. The temperature T aero called aerodynamic temperature represents the average temperature of the air in the vicinity of the vegetation elements within the canopy, at the height of the aerodynamic level (defined as the sum of the displacement height and the roughness length for momentum). There is no measuring device for this temperature, which is usually inverted from turbulent flux measurements. Moreover, it can be significantly different from the ensemble directional radiometric surface temperature T surf (as defined by [START_REF] Norman | Terminology in thermal infrared remote-sensing of natural surfaces[END_REF] which is usually derived from brightness temperature measurements made by a thermo-radiometer or infrared thermometer at nadir or at a specific view angle [START_REF] Kustas | Advances in thermal infrared remote sensing for land surface modeling[END_REF][START_REF] Kustas | Utility of radiometricaerodynamic temperature relations for heat flux estimation[END_REF][START_REF] Stewart | Sensible heat-flux radiometric surface-temperature relationship for 8 semiarid areas[END_REF] and which is used to assess energy balance from remote sensing data. In general, the relationship between aerodynamic and surface temperatures is obtained through the use of a dual-source energy balance when either the vegetation or the soil bulk skin temperatures are known [START_REF] Lhomme | Sensible heat flux-radiometric surface temperature relationship over sparse vegetation: Parameterizing B-1[END_REF].

Retrieving surface temperature for each of the different components of the surface (mostly soil and plants) is difficult, especially with current remote sensing platforms [START_REF] Jia | A practical algorithm to infer soil and foliage component temperatures from bi-angular ATSR-2 data[END_REF].

The use of single-source models is therefore favoured over dual-source models to estimate pixel average turbulent fluxes from a mixed-pixel radiometric temperature. For those models, there is a need to develop robust yet simple methods to relate the aerodynamic temperature to the surface temperature. This has been subject of debate for a long time (e.g. see Carlson et n.4 al., 1995;[START_REF] Kalma | Estimating evaporation from pasture using infrared thermometry: evaluation of a one-layer resistance model[END_REF]. Many formulations exist in the literature, and a comprehensive terminology and conversion formulae are proposed by [START_REF] Matsushima | Relations between aerodynamic parameters of heat transfer and thermal-infrared thermometry in the bulk surface formulation[END_REF].

Historically, most expressions governing the relationship between the aerodynamic and the surface temperatures have been built on an analogy between wind and temperature profiles within the canopy. However, while the bottom boundary condition for wind (a null value at the height of the aerodynamic level) allows defining a roughness length for momentum exchange (z om ), the bottom value of the temperature profile, the aerodynamic temperature, is generally unknown. One assumes usually that a roughness length z or (improperly named roughness length for thermal exchange) can be defined so that, at a height corresponding to the displacement height plus z or , the surface temperature of the vegetation can be considered as representative of the aerodynamic temperature. The relationship between both roughness lengths translates into what Matsushima names the "radiometric kB -1 " (kB -1 radio ) which is written as

1   radio kB om or e z z
. In that case, the difference between the aerodynamic temperature and the surface temperature in the sensible heat flux formulation is expressed by an additional resistance [START_REF] Lhomme | Sensible heat flux and radiometric surface temperature over sparse Sahelian vegetation .2. Model for the kB(-1) parameter[END_REF]:

                                     L z L d z z d z T T ku c H ku kB r T T c H
                                
is the aerodynamic resistance for heat exchange before kB -1 correction, u * the friction velocity, k the von Karman constant, z the measurement height of the atmospheric forcing, d the displacement height, L the Monin Obhukhov length and  h the stability correction function for heat transfer given by [START_REF] Paulson | The Mathematical Representation of Wind Speed and Temperature Profiles in the Unstable Atmospheric Surface Layer[END_REF]. kB -1 radio is derived according to the expected air temperature profile within the canopy and expressed as a function of meteorological data, LAI and plant height. Amongst the well n.5

known formulae, one can cite those from Blümel [START_REF] Blümel | Estimation of sensible heat flux from surface temperature wave and onetime-of-day air temperature observation[END_REF], Massman [START_REF] Massman | A model study of kB(H)(-1) for vegetated surfaces using 'localized near-field' Lagrangian theory[END_REF] revisited by [START_REF] Su | An evaluation of two models for estimation of the roughness height for heat transfer between the land surface and the atmosphere[END_REF] and Lhomme [START_REF] Lhomme | Sensible heat flux-radiometric surface temperature relationship over sparse vegetation: Parameterizing B-1[END_REF].

Other authors have proposed a somewhat simpler, and easier to interpret, formulation of the relationship between T surf and T aero , called the " function", originally proposed by [START_REF] Chehbouni | Estimation of sensible heat flux over sparsely vegetated surfaces[END_REF].  is expressed solely in terms of the temperatures, independently from wind speed: Even for isothermal surfaces, usually bare soils or very dense canopies, the aerodynamic temperature can be slightly different from the surface temperature, because the diffusion process for heat transfer adds to the convective exchange of air. There is therefore a difference between the effective eddy diffusivities for momentum and heat exchange, which can be again translated into an excess resistance function of an "aerodynamic kB -1 " or kB -1 aero . The available kB -1 formulae, derived either empirically or from scalar and flux theoretical profiles in the canopy, account for both aspects: the difference between the aerodynamic and the radiometric temperatures (radiometric kB -1 ) and the difference between momentum and heat exchange diffusion processes (aerodynamic kB -1 ). For isothermal surfaces, one can expect that there is no radiometric component within the combined (radiometric and aerodynamic) kB -1 . The combined kB -1 retrieved for those surfaces from observations by solving for kB -1 radio in Eq. 3 is fairly low (within the range 0-5 according to [START_REF] Verhoef | Some Practical Notes on the Parameter kB -1 for Sparse Vegetation[END_REF][START_REF] Massman | A model study of kB(H)(-1) for vegetated surfaces using 'localized near-field' Lagrangian theory[END_REF][START_REF] Yang | Turbulent flux transfer over bare-soil surfaces: Characteristics and parameterization[END_REF]. For strongly non-isothermal situations, which correspond in general to n.6

h a air surf p r T T c H ,     (4) i.e.
intermediate LAI values (LAI in the range 0.5-2), the combined kB -1 is much higher (in the range 10-30 according to the same authors). One can thus assume that kB -1 aero is usually smaller than kB -1 radio for all LAI values, or that the difference between the surface temperature and the aerodynamic temperature will have on average a much larger impact on the sensible heat flux than the difference between the diffusion processes for heat and momentum at the vicinity of the canopy. While the radiometric kB -1 do not discriminate between both differences (the difference between the aerodynamic and surface temperatures and the difference between the roughness lengths for momentum z om and heat exchange z oh ), the  function allows us to separate both aspects and keep the difference between z oh and z om in the formulation of the aerodynamic resistance : Consequently, Eq. 4 can be rewritten as:

                                     L z L d z z d z T T ku c H ku kB r T T c H oh h h oh air surf p aero hm a air surf p ln 1 ,     (6)
One must also note that the published values of the radiometric kB -1 (kB -1 radio ) [START_REF] Matsushima | Relations between aerodynamic parameters of heat transfer and thermal-infrared thermometry in the bulk surface formulation[END_REF] according to Lhomme, Blümel and Massman/Su (see formulations in Table 1) can be converted into  values by combining equations 3 and 6:

       ku kB r ku kB r radio hm a aero hm a 1 , 1 ,  (7)
The Blümel and Massman/Su formulations depend on Leaf Area Index LAI through, mostly, the fraction cover (f c ) of the canopy and on two parameters difficult to assess, the component aerodynamic kB -1 for bare soil and for vegetation canopy, respectively (See Table 1). Since for bare soil and full cover conditions there is no large difference between both temperatures, n.7 the  function is fairly easily interpreted:  values are close to 1 for those bare soil and full cover conditions, that is, more generally, for all homogeneous isothermal surfaces, while for sparse vegetation  decreases. In those conditions, the soil temperature has a large impact on the radiative surface temperature whereas the aerodynamic temperature remains closer to a mix of air and vegetation temperatures and is less influenced by the soil temperature. Since the soil temperature around midday is generally higher than the vegetation temperature, the observed radiative surface temperature is often larger than the aerodynamic temperature around that time. Factors influencing  and kB -1 radio include LAI and other plant geometrical features such as height and fraction cover, friction velocity, time of the day, solar radiation etc.

However, most studies agree on the fact that LAI is by far the main driving factor, at least for agricultural canopies for which the turbid medium (random leaf dispersion, [START_REF] Myneni | A review on the theory of photon transport in leaf canopies[END_REF] and permeable-rough transfer hypotheses are valid [START_REF] Kustas | Utility of radiometricaerodynamic temperature relations for heat flux estimation[END_REF][START_REF] Verhoef | Some Practical Notes on the Parameter kB -1 for Sparse Vegetation[END_REF]. This is further confirmed by dual-source land surface models which predict the aerodynamic temperature through the classical dual-source approach [START_REF] Shuttleworth | Evaporation from sparse crops -an energy combination theory[END_REF].

The evolution of as a function of LAI is presented in Figure 1a as obtained from an uncalibrated run of the ICARE SVAT model [START_REF] Gentine | Analysis of evaporative fraction diurnal behaviour[END_REF] for the B124 wheat site in Morocco (see below). One can observe that the simulated shape of the (LAI) relationship decreases sharply from 1 when LAI increase to about 1, and increases more slowly for higher LAI, tending again to 1 for LAI well above 3. The lognormal distribution function is therefore a good candidate to represent the evolution of 1-(LAI) for the whole range of LAI values.

Consequently, we propose the following empirical relationship for  :

  2 2 * 2 ) ln( 2 * * 1 b c LAI e b LAI a        (8)
where a, b and c are empirical coefficients that need to be calibrated.

n.8

The objectives of the present paper are threefold:

1-to retrieve variations with LAI from observations for nine experimental cultural cycles where seasonal evolution of factors governing  is available and LAI values range from 0 to well above 2, and by doing so, assess the variability in shapes and scales of the(LAI) relationship, 2-to compare several formulae of (LAI), including the new one (Eq. 8; hereafter referred to as the Boulet et al. expression) , against observed trends and 3-to compare the performances of the various formulae in computing sensible heat flux from observed in-situ radiometric surface temperature.

The new Boulet et al. expression (Eq. 8) of the  function will be first calibrated over the values of  derived from experimental data acquired on the B124 wheat site in Morocco (i.e. values will be obtained for a, b and c in Eq. 8). Next, it will be tested on data acquired over eight other crop cycles in South of France.

STUDY SITES

All study sites cover an entire agricultural season, from bare soil to harvest.

The first dataset was collected at the B124 site (31.67250°N, 7.59597°W) in the R3 irrigation perimeter in the Haouz semi-arid plain in Morocco during the SudMed project [START_REF] Chehbouni | An integrated modelling and remote sensing approach for hydrological study in arid and semi-arid regions: the SUDMED programme[END_REF][START_REF] Duchemin | Monitoring wheat phenology and irrigation in Central Morocco: On the use of relationships between evapotranspiration, crops coefficients, leaf area index and remotely-sensed vegetation indices[END_REF] in 2004. The climate is semi-arid with an average annual rainfall of the order of 150 mm. The chosen field (number B124) was cultivated with winter wheat and its size (4 ha) exceeded the basic fetch requirements. LAI and vegetation height ranged from 0 in January to 4.5 and 0.8m at maximum development (April), respectively (See [START_REF] Boulet | Monitoring water stress using time series of observed to unstressed surface temperature difference[END_REF], for more information on this dataset).

n.9

The second dataset was collected over two cultivated plots, Auradé (43°54'97''N, 01°10'61''E) and Lamasquère (43°49'65''N, 01°23'79''E), separated by 12 km and located near Toulouse (South West France). The climate is Mediterranean with an average annual rainfall of the order of 620 mm for both 2006 and 2007, which can be considered as "average" years as far as rainfall is concerned. The Auradé plot was cultivated with winter wheat Irrigation was applied in particular to the Sorghum and Peas crops.

n.10

All sites were equipped with a tower where standard meteorological forcing data were acquired. Plant height and total (green and dry) Leaf Area Index were measured every month or so by planimetry and hemispherical photography and the resulting values have been interpolated to daily estimates. Half hourly sensible heat fluxes were measured with a Campbell Sci. CSAT3 (B124, Auradé and Lamasquère sites) or a Young 81000 (Avignon site) 3D sonic anemometer. Auradé, Lamasquère and Avignon are part of the CarboEurope-IP Regional Experiment [START_REF] Dolman | The CarboEurope regional experiment strategy[END_REF] and the CarboEurope-IP Ecosystem Component.

In that context, the data were used for analyzing CO 2 surfaceatmosphere exchanges and production for fields with a large interannual rotation of crop types (e.g. [START_REF] Kutsch | The net biome production of full crop rotations in Europe[END_REF].

For those sites, the Level 3 flux products (i.e. non gapfilled) were used.

For all sites, TIR data were acquired with a nadir looking 60° Field Of View Apogee IRTS-P (R3 and SudOuest sites) or an Heitronics KT15 (Avignon site) Infra Red Thermoradiometer and a Kipp and Zonen CNR1 hemispherical pyrgeometer. For the R3 B124 site, the IRTS-P device has been calibrated using an Everest black body during the experiment and prior to the experiment in a laboratory with an adjustable ambient temperature. Retrieved surface temperatures from both instruments showed a bias of 0.21°C and a Root Mean Square Difference of 1.16°C for instantaneous values at midday. Surface temperature estimates from the CNR1 were used for this site. The KT15 instrument in Avignon was calibrated in the laboratory calibration facilities by looking at a black-body at various temperatures. An error analysis of the whole calibration-measurement chain gave an accuracy of 0.4°C.

PROCEDURE FOR RETRIEVING 

                                         ( 10
)
Where L obs is the Monin Obhukhov length estimated from observations and u *,obs the measured friction velocity.

In order to assess either  (from Eq. 8) or the radiometric kB -1 , one would need in theory a full determination of the three main variables: z om , d and z oh =z om /exp(kB -1 aero ). This is feasible from single-direction radiometric measurements only for bare soil and full cover conditions since in the latter cases aerodynamic and surface temperatures can be considered as identical (kB -1 radio =kB -1 aero or =1). For all other conditions, z om and z oh can be retrieved from u * and scale temperature

* * u c H T p  
measurements only if T aero is already well known, i.e. if the radiometric kB -1 is known. As a consequence, we decided to use the following expressions from [START_REF] Colaizzi | Comparison of aerodynamic and radiometric surface temperature using precision weighing lysimeters[END_REF] and [START_REF] Pereira | Evapotranspiration: Concepts and future trends[END_REF] to compute z om and d:

            LAI e h d LAI 5 . 0 1 2 1 (11)   LAI LAI om e he z 5 . 0 5 . 0 1     (12)
where h is the vegetation height.

n.12

Given d, one can derive an estimated roughness length z om,obs according to:

                 obs m obs obs L d z u ku obs om e d z z , , (13)
where u obs is the wind velocity at reference height. We checked that z om values obtained by Eq. 12 are valid, i.e. consistent with z om,obs obtained by Eq. 13 (not shown).

We choose emissivity values (between 0.96 for a bare soil and 0.98 for a vegetation at full cover) for the computation of surface temperature from brightness temperature measurements and a null value for the aerodynamic kB -1 (

1 1   aero kB e
) so that retrieved  obs values deduced from Eq. 9 tend to one for very high LAI values (and therefore retrieved radiometric kB -1 will be close to zero when solving for kB -1 radio in Eq. 7), which is what is expected for full cover conditions. Values close to zero for the radiometric kB -1 can be found elsewhere in the literature, notably for quasi-bare soils, as suggested both by [START_REF] Massman | A model study of kB(H)(-1) for vegetated surfaces using 'localized near-field' Lagrangian theory[END_REF](in his paper the radiometric kB -1 tends to 0 when LAI=0 or LAI>>2) and [START_REF] Matsushima | Relations between aerodynamic parameters of heat transfer and thermal-infrared thermometry in the bulk surface formulation[END_REF] in spite of the large Reynolds numbers encountered in typical agricultural bare soil situations. Of course, overall retrieval results and performances depend strongly on the accuracy of the surface temperature measurements, and this should be investigated further, but it is beyond the scope of this short paper. On the other hand, neglecting kB -1 aero also leads to the best performance in estimating sensible heat flux for all LAI values and all sites (not shown) which justifies our choice. [START_REF] Yang | Turbulent flux transfer over bare-soil surfaces: Characteristics and parameterization[END_REF] found values of the aerodynamic kB -1 for bare soils between -5 (at night) and 5 (around midday) with an average between 0 and 5 for positive values of observed sensible heat flux. In what follows,  obs values are computed for unstable conditions selected on the basis of the following rules: H obs >0, total incoming solar radiation >10 W/m 2 and T surf,obs >T air,obs . Retrieved  obs values from Eq. 9 between -1 and 2 are kept since they represent realistic values of the temperature profile within the canopy and the soil/vegetation n.13 mixed/ensemble contribution (this represents >85% of all H obs data between 10am and 4pm

for the R3 B124 dataset). Blümel, Massman/Su and Lhomme expressions, when LAI tends towards 0, the turbulent behavior gives more weight to the aerodynamic kB -1 for a bare soil according to [START_REF] Brutsaert | Evaporation into the atmosphere[END_REF] which tends to underestimate  [START_REF] Yang | Turbulent flux transfer over bare-soil surfaces: Characteristics and parameterization[END_REF]. As a consequence, the sensible heat flux simulated from the different formulae for using the observed surface temperature (Equation 4), as one would do with remote sensing observations, is closer to the observed sensible heat for the proposed new formula compared to the three previous ones, especially for quasi bare soil (LAI<0.5) conditions. This appears clearly on Figure 2, where the resulting sensible heat flux is estimated by using the various formulae for kB -1 or  together with Eqs. For all sites, the new formulation (Eq. 8) fits fairly well the observed , except for Lamasquère in 2007, which means that the calibrated values for a, b and c at the R3 B124 site are well suited for other sites. For sunflower, corn and sorghum however (Figure 3c, e and f), the trend in (LAI) matches the observed trend but not the observed amplitude, and the value for a should therefore be lower for those cover types. This might be due to the size and shape of the leaves, or to the fact that these canopies show more defined geometrical features (rows, preferred orientation of leaves and flowers…) and are less easily described by turbid medium (for radiation) or permeable-rough interfacial layer (for turbulent transfer) theories. It is also possible to adjust the a, b and c values to fit more closely the observed (LAI) curve. Adjusted values for a range between 1.2 and 1.7, and adjusted values for b and c between 0.5 and 0.8, but do not translate into a significant improvement in computing sensible heat flux.

Since the primary objective of the  formulation is to provide accurate estimates of sensible heat fluxes from observed surface temperature, one needs to assess the resulting performance in using the various  models to simulate H. However, it should be noted that the accuracy of those estimates depends primarily on the precision on surface and air temperatures, wind speed, roughness length and the validity of the Monin-Obukhov Similarity Theory.

The performance in estimating sensible heat flux for all sites was assessed using the various  formulations. Results are shown in Table 2. All formulations perform well for some sites and much less well for others, but the new formulation shows comparable or better performances than the others as indicated by the small number of RMSE values above 60 W/m 2 as well as the high number of values under 50 W/m 2 . The simplicity of this formulation and the fact that it allows tuning at least one parameter (a) to provide realistic values of the difference between the aerodynamic and the surface temperature are two advantages of this formulation. One n.16 must note that, for both the Blümel and Massman formulations, the poorly known parameters corresponding to the aerodynamic kB -1 for bare soil and full cover (which are found to be close to zero here) can be adjusted as well as the parameter of the exponential decay of the Beer Lambert law used to convert LAI to fraction cover, but in our simulations the results did not prove to be very sensitive to the latter (not shown). Again, all formulations perform well for canopies well described by the turbid medium theory, and less for plants with more defined geometrical features or larger intercropping patterns [START_REF] Barillot | How good is the turbid medium-based approach for accounting for light partitioning in contrasted grass-legume intercropping systems ?[END_REF].

CONCLUSION

An empirical formulation of the difference between the aerodynamic and surface temperatures as a function of Leaf Area Index has been proposed, which represents in a realistic way the observed variations and leads to satisfactory performance in simulating the sensible heat flux compared to other existing formulations. It should be noted though that the observed variations in this difference were assessed using a null aerodynamic kB -1 , based on the fact that the radiometric kB -1 should be close to zero (or  close to one) for bare soil and full cover conditions. However, the assumption that the difference between the surface temperature and the aerodynamic temperature will have on average a much larger impact on the sensible heat than the difference between the diffusion processes for heat and momentum at the vicinity of the canopy should be investigated more thoroughly since there is no agreement on what value should be used for the aerodynamic kB -1 for bare soil, intermediate cover and fully covering vegetation, respectively [START_REF] Verhoef | Some Practical Notes on the Parameter kB -1 for Sparse Vegetation[END_REF]. In particular, there is no evidence whatsoever that the roughness length for heat exchange should be a constant fraction of the roughness length for momentum for all LAI values, as it is commonly assumed in SVAT models.

n.17

(

  Triticum aestivum L., maximum LAI: 3.8, maximum height: 0.68 m) from Oct-2005 to Jun-2006 and with sunflower (Helianthus annuus L., maximum LAI: 2.0, maximum height: 1.27 m) from Apr-2007 to Sept-2007. The Lamasquère plot was cultivated with maize (Zea mays L., maximum LAI: 3.3, maximum height: 2.3 m) used for silaging from May-2006 to Aug-2006 and with winter wheat from Oct-2006 to Jul-2007 (maximum LAI: 4.5, maximum height: 0.83 m). This site was irrigated in 2006 when maize was cultivated. For a complete description of the site characteristics, management practices, biomass inventories, vegetation area measurements, instrumentation setup and fluxes calculation procedures see Beziat et al. (2009). The third dataset was acquired on the Avignon 'Flux and Remote Sensing Observation Site', located in Provence, in southeastern France (N 43,92°; E 4,88°; altitude 32 m). The region is also characterized by a typical Mediterranean climate (annual climatic mean of 14° C for temperature and of 680 mm for precipitation). However, during the 2004-2007 period the yearly average for rain was 450 mm with a high variability (313 mm -745 mm). During the observational period the crop rotation was: Durum Wheat (Triticum durum, maximum LAI: 4.0, maximum height: 1.0 m) from January till June 2004, Peas (Pisum sativum, maximum LAI: 2.9, maximum height: 0.43 m) from April till June 2005, Durum Wheat (maximum LAI: 5.5, maximum height: 0.75 m) again from November 2005 to June 2006 and Sorghum (Sorghum bicolor, maximum LAI: 3.0, maximum height: 1.16 m) from May to August 2007.

For

  all sites, retrieved  function values estimated from observations,  obs , are derived from measured sensible heat H obs , surface temperature T surf,obs and meteorological data such as air n.11 temperature T air,obs that influence the aerodynamic resistance r a,obs calculated using the Moninr a,obs is computed using measured friction velocity values. Therefore:

  4. RESULTS FOR WHEAT AND OTHER CROPS 4.1. Results for the R3 B124 wheat site (calibration of Eq. 8) Median values of the scattered  obs values are shown for each 0.5 LAI interval on Figure 1a, together with  values simulated by the ICARE dual source SVAT model applied to the R3 B124 site without any calibration of ICARE,  values interpolated along LAI values from individual  estimates proposed by Matsushima (2005) as well as values obtained from Eq. 8 (also referred to as "Boulet et al."), whose a, b and c coefficients are manually adjusted to fit  obs values: a=1.7, b=c=0.8. Error bars are shown for observed  values by assuming for Eq. 9 an error of 1K on T surf,obs , 10% on u *,obs and z om and 10 W/m 2 on H obs . values converted from Eq. 7 based on the three kB -1 expressions proposed by Blümel, Massman/Su and Lhomme are also shown in Figure1a.For well developed vegetation (say, above LAI=2), all expressions are fairly similar or at least show consistent trends. For sparse and less developed vegetation however (LAI<1, figure1b, where all single "observed" retrieved  values are shown instead of median values), there is a large difference between the Blümel, Massman/Su, Matsushima and Lhomme expressions, on the one hand, and ICARE outputs, the proposed Boulet et al. expression and the observations, on the other hand. Note that ICARE has not been calibrated to the R3 dataset, therefore only the overall shape of the (LAI) relationship is being used to introduce Eq. 8, not its particular scaling/extent properties (a, b and c values) which are derived from  obs values. The differences between both groups of estimates could be explained by the fact that for the n.14

  3 and 4 respectively are shown. On this figure, the sensible heat flux values away from the 1:1 line for the Blümel, Massman/Su, Matsushima and Lhomme expressions correspond to those conditions (LAI<0.5).One must note that if we ignore the  correction (i.e. by setting =1 or root mean square error RMSE between observed and simulated sensible heat fluxes increases from ~50 to ~80 W/m 2 for the entire season. This confirms that taking into account the difference between the aerodynamic and the surface temperature is crucial to deriving accurate turbulent heat fluxes.4.2. Results for 8 agricultural seasons in Southern France (validation).The different expressions for  are tested for the 3 other sites, i.e. 8 growing seasons and different vegetation and climate conditions. We assume that the adjusted values of the three parameters for the new formulation (a=1.7, b=c=0.8) can provide a good estimate of (LAI) fluctuations for agricultural areas. These values are therefore kept for all sites. Results are shown in Figure3. Surface temperature is estimated either from the hemispherical (CNR1) or the directional (KT15, IRTS-P) sensors, as specified in the captions. In general, the hemispherical device produces a (LAI) relationship closer to the expected concave-up shape.

ACKNOWLEDGEMENTS

Fundings from the Centre National d'Etudes Spatiales (CNES) for the MiSTIGrI (MicroSatellite for Thermal Infrared Ground surface Imaging) phase A study, the MISTRALS (Mediterranean Integrated STudies at Regional And Local Scales) SICMed (Continental Surfaces and Interfaces in the Mediterranean area) program, and the SIRRIMED (Sustainable use of irrigation water in the Mediterranean Region) European project (FP7 -grant agreement 245159), are gratefully acknowledged. Em. Prof. Jetse D. Kalma is also thanked for his very insightful comments on an earlier draft of the manuscript. Avignon and SudOuest data were acquired and processed in the frame of the CARBOEUROPE-IP and the CARBOFRANCE project funded by the European FP7 Program (GOCECT-2003-505572) and the French Ministry in charge of Environment (GICC programme). For the R3 dataset, in addition to IRD, financial support was provided by EC in the frame of the WATERMED project (contract ICA3-CT-1999-00015) and IRRIMED project (contact ICA3-2002-10080) and by the French Programme National de Télédétection Spatiale (PNTS) and the French space agency (CNES).