domain: spi phys

text:
Diogo Sabino
email: diogo.ferreirasabino@onera.fr

Olivier Marquet

David Fabre

Vincent Mons

Global Aeroelastic Stability Analysis of a NACA0012 Airfoil in Transitional Reynolds Regime

Keywords: linear stability analysis, transitional Reynolds regime, laminar separation flutter

The onset of pitch-oscillations of a NACA0012 airfoil mounted on a torsion-spring is investigated by the means of a global fluid-structure linear stability analysis (LSA). Experimental studies from literature have shown this configuration undergoing on a flutter behaviour for a close range of Reynolds. This phenomenon is attributed to the laminar character of the boundary layer on the airfoil surface and to its detachment, transition and subsequent reattachment. In this regard, we carried out a global fluid-structure LSA of a mean flow field, issued from a time and spanwise averaged Direct Numerical Simulation, at different incidences, α, for Re = 50 000. The LSA at α = 0 • predicts a static divergence behaviour of the airfoil from its equilibrium position, whereas the LSA at different incidences predicts the stabilisation of the static mode, associated to the existence of a laminar separation bubble over the upper surface of the airfoil.

Introduction

The linear stability of fluid-solid systems is a classic subject on the aeroelasticity community [START_REF] Bisplinghoff | Aeroelasticity[END_REF]. The primordial aim was to delimit the flight envelope of military and commercial aircraft cruise flights, whose Reynolds number often exceeds a million and where the boundary layer remains turbulent and attached to the airfoil for a relatively hight range of angles of attack. However, as the Reynolds number decreases, the transition location moves downstream, originating a laminar boundary layer over a considerable portion of the wing. For a low-to-moderate Reynolds number range, 10 4 ≤ Re ≤ 10 6 , the flow regime is characterised by the coexistence of laminar, transition and turbulent regions, where laminar boundary layer region can lead to complex viscous phenomena, such as laminar boundary layer separation, with a possible transition and reattachment. In this Reynolds regime, Poirel et al. [START_REF] Poirel | Self-sustained aeroelastic oscillations of a NACA0012 airfoil at low-to-moderate Reynolds numbers[END_REF] carried out a wind tunnel experimental investigation of a NACA0012 airfoil mounted on a torsional spring system, demonstrating the existence of a self-sustained limit cycle oscillation (LCO) within the Reynolds number range of 45 × 10 3 ≤ Re ≤ 130 × 10 3 . This limit cycle was characterised by a Strouhal frequency around St = 0.06 and an oscillation amplitude of α max ≈ 5.5°, around the equilibrium position at α = 0°. The authors suggested that the loss of stability of the airfoil at α = 0°and, consequently, the dynamics of the pitch oscillations, were governed by the laminar separation of the boundary layer, either due to the trailing edge separation or to the presence of the laminar separation bubble (LSB). The phenomenon was numerically reproduced by [START_REF] Poirel | Computational aeroelastic simulations of self-sustained pitch oscillations of a NACA0012 at transitional Reynolds numbers[END_REF] and others, using a time-marching coupled aeroelastic simulation, where it was confirmed that the laminar separation near the trailing edge plays a critical role in initiating and sustaining the pitch oscillations. For this reason, the phenomenon was labelled as Laminar Separation Flutter (LSF). Regarding the global linear stability analysis (LSA), the aim is to evaluate the linear stability character of the dynamical system considered. In the case of coupled fluid-structure systems, several configurations have been studied in the past years [START_REF] Assemat | The onset of unsteadiness of two-dimensional bodies falling or rising freely in a viscous fluid: a linear study[END_REF][START_REF] Sabino | Vortex-induced vibration prediction via an impedance criterion[END_REF]. Concerning the LSF, Negi [START_REF] Negi | Stability and Transition in Pitching Wings[END_REF] performed a fluid-structure LSA over a NACA0012 airfoil at Re = 50 000 and α = 0°around a mean flow issued from a spanwise and time-averaged time-marching simulation. They found an unstable static mode, responsible for the departure of the structure from its equilibrium position, and further associate the subsequent non-linear oscillation to a sub-critical mode, which would be non-linearly excited. In this work, we propose to give an insight on the linear behaviour of the mean flow field at Re = 50 000 by the means of a global fluid-structure linear stability analysis for different incidences. In a first time, the analysis is carried out at α = 0 • . This analysis is then followed by a discussion on the mean flow topology as the angle of attack increases, and in a second time, we analyse the cases for α > 0°.

Mathematical and Numerical Formulations

Fluid and Solid Models

The dynamics of an incompressible flow field is governed by the Navier-Stokes equations, that can be posed, with respect to a non-inertial frame of reference, as

∀X ∈ Ω f ⊂ R 3 , t > 0 :      ∂U ∂t +Ω × U +([U -u w ] •∇) U + ∇p - 1 Re ∇• ∇U + (∇U ) T = 0 , ∇•U = 0 , (1a) (1b)
with

Ω f the fluid domain, U (X, t) = [U, V, W ]
T the absolute fluid velocity field, p(X, t) the fluid pressure and u w the airfoil velocity, defined as

u w = Ω × (X -X EC ) + ẊEC with ẊEC = U ∞ [-cos(θ), sin(θ), 0] T . (2)
The Reynolds number is defined as Re = cU ∞ /ν, with c the airfoil chord, U ∞ its the translation velocity on the absolute frame and ν the kinematic viscosity. The airfoil is considered as a rigid body, allowed to rotate around its elastic centre, X EC , fixed at 0.186c. The pitching motion is characterised by the instantaneous angular displacement, θ(t), and by the angular velocity, Ω(t), both considered positive for a nose-down displacement. Ω(t) represents the e z -component of the vectorial angular velocity Ω = Ω e z . Additionally, the angle of attack is defined simply as α = -θ. We consider homogeneous Dirichlet conditions on the inlet and lateral domain boundaries for the velocity field and a natural Neumann condition at the outlet, imposing a stress-free condition.

The fluid state vector can be posed as

q f = [U , p] T .
The Direct Numerical Simulations (DNS) were performed for a fixed airfoil configuration, hence having Ω = 0 and θ fixed on time. The motion of the airfoil is considered for the linear stability analysis formulation, with a temporal evolution defined as

I s (X EC ) dΩ dt + D s Ω + K s θ = m z q f , X EC , with dθ dt -Ω = 0 (3)
where

I

Linearised Equations

The linear stability of the airfoil structure is determined by investigating the nature of small amplitude perturbations of a steady state field. In that sense, the the fluid-solid state variable, q, is decomposed as the sum of a steady solution q = q f , θ, Ω T and unsteady infinitesimal perturbation q as q(X, t) = q(X) + q (X, t) ,

where 1. As the steady solution is independent of time, its angular velocity is equal to zero, Ω ≡ 0, unlike the angle of incidence θ = 0, which can represent the airfoil at a given equilibrium position. The infinitesimal perturbation can be further decomposed into the form of global modes, i.e.,

q (X, t) = q(X) e σt + c.c. 2 , (5)
where c.c. is a short hand for complex conjugate, q(X) is a complex eigenmode with an eigenvalue σ = λ + iω. The eigenmode describes the spatial distribution of the perturbation, whereas the eigenvalue accounts for amplification rate, λ, and frequency characteristics, ω, for that eigenmode. The sign of the amplification rate indicates if the corresponding eigenmode decays (negative sign) or grows (positive sign) exponentially in time, dictating if the steady state is linearly stable or unstable to the eigenmode perturbation. The complex moment of the eigenmode is defined as mz = m T qf . The steady state is here represented by the mean flow field, issued from the time and spanwise averaged DNS, at different incidences, for Re = 50 000. As the mean flow is averaged in the spanwise direction, the LSA computations are made in a two-dimensional space, having q = Û , V , p T .

Injecting the decomposition 4 into the equations 1 and 3, one obtains, at order , an eigenvalue problem, that can be put in the form σMq = Lq:

σ     M f 0 0 0 1 0 0 0 I s         qf θ Ω     =     L f v θ v Ω 0 0 1 m T -K s -D s         qf θ Ω     , (6)
where M f is the fluid mass operator, L f the fluid Jacobian operator, and v θ , v Ω the vector components of the solid motion. v θ represents an unitary linear displacement of the airfoil structure, whereas v Ω represents an unitary linear angular velocity.

Numerical Formulation

The spatial discretisation is based on a finite element method. The unknown velocity and pressure fields were discretised using the classical Taylor-Hood basis in order to satisfy the Ladyzhenskaya-Babuška-Brezzi (known as LBB) condition (see, for instance, [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF]). All the discrete matrices resulting from the projection of the variational formulations onto the basis of finite elements were built with the FreeFem++ software [START_REF] Hecht | New development in freefem++[END_REF], and the linear systems solved via a parallel implementation using SLEPc/PETSc software [START_REF] Hernandez | SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems[END_REF][START_REF] Balay | PETSc users manual[END_REF]. The resolution of the shift-and-invert system passed through an Arnoldi iterative algorithm [START_REF] Saad | Numerical Methods for Large Eigenvalue Problems[END_REF].

3 Results

Fluid Stability Analysis at

Re = 50 000, α = 0 •
We start the analysis by a purely fluid eigenproblem, at Re = 50 000, where the airfoil is fixed at α = 0°. The eigenvalue results are presented in the one-sided spectrum in figure 1a, with the growth rate on the horizontal axis and the eigenmode frequency on the vertical axis. A zoom at the origin is plotted in figure 1b. The unstable zone, in grey, corresponds to a positive amplification rate. The spectrum is characterised by a cloud of stable modes and a single unstable eigenmode, vibrating at a frequency of ω = 29.73. This frequency corresponds to the frequency of the vortex-shedding that can be found on the unsteady airfoil wake. Although this mode is expected to be marginally stable, as in the case of the unsteady vortex shedding on a cylinder wake at supercritical Reynolds regimes [START_REF] Sipp | Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows[END_REF], the discrepancy can be associated to the high convective character of the perturbation. This fact is in accordance with the conclusions of [START_REF] Poirel | Self-sustained aeroelastic oscillations of a NACA0012 airfoil at low-to-moderate Reynolds numbers[END_REF], when they found no correlation between the pitch response frequencies and the vortex-shedding frequency and in accordance with [START_REF] Negi | Stability and Transition in Pitching Wings[END_REF], where the projection of the solid variables onto the fluid variables was concluded to be close to zero for all the high frequency eigenmodes.

FSI Stability Analysis at

Re = 50 000, α = 0 •
The FSI stability analysis was carried out for the same dimensional parameters as the ones reported in [START_REF] Poirel | Self-sustained aeroelastic oscillations of a NACA0012 airfoil at low-to-moderate Reynolds numbers[END_REF] and subsequently used by [START_REF] Negi | Stability and Transition in Pitching Wings[END_REF], i.e., I dim s = 0.00135kg m 2 , D dim s = 0.002N m s and K dim s = 0.15N m/rad. Figure 1c shows the one-sided spectra of the FSI linear stability analysis. Differently from the hydrodynamic case, the FSI eigenproblem reveals the emergence of two pairs of eigenvalues, associated to the structure. The first pair is located along the real axis, presenting no associated frequency. The unstable eigenmode has a growth rate of λ = 0.2091, whereas the other eigenvalue is stable and located in the vicinity of the fluid modes. The second pair is located at the a growth rate of λ = -0.1500 and with a frequency of ω = ±0.5264. The eigenmode presenting a divergent character will be labelled Static Mode, whereas the low frequency eigenmode, possibly at the origin of the nonlinear flutter oscillations, will be labelled Flutter Mode. The flutter mode of the fully coupled fluid-structure problem can be compared with the eigenvalues of the structure equation. The real and imaginary part of these eigenvalues can be computed analytically and read

σ 2 s I s + σ s D s + K s = 0 ⇒ λ s = - D s 2I s ω s = ± K s I s -λ 2 s . (7)
From this result, one can infer that a structure without damping will be marginally stable, resonating at the natural frequency of ω s = (K s /I s ). Any positive value of the damping coefficient results into a negative growth rate, damping any structure oscillation, which returns to a stable equilibrium position. For the current parameters, the natural structure eigenvalues are σ s = -0.0234 ± 0.3320i. The positive frequency eigenvalue can be visualised in the one-sided spectrum of figure 1c and compared to the FSI flutter eigenvalue. The coupled fluid-structure problem gives a more stable mode, suggesting that, for this incidence, the fluid has a stabilising effect on the structure equation.

FSI Stability Analysis at

Re = 50 000, α > 0 •
The DNS computations were further carried out at the incidences from α = 0.5 • to α = 2.0 • , in order to access to the time and spanwise averaged mean flow field at angles of attack higher than 0 • . The streamwise component of the velocity field of the mean flow for four angles of attack is present in the figure 2, at left. The upper figure corresponds to the case at α = 0 • . As illustrated by the streamlines, the flow field has a separation point at x/c = 0.7. This separated flow gives rise to a long recirculation bubble that is present in the rear of airfoil, in both lower and upper surfaces. As the angle of attack increases, the separation point on the upper surface of the airfoil moves upstream, while the one on the lower surface moves downstream, the former

0.0 0.2 -0.2 0.0 0.4 0.8 λ ω α = 0.0 • α = 0.5 • α = 1.0 • α = 1.5 • (a) 0.0 0.2 -0.2 0.0 0.4 0.8 λ ω α = 0.5 • α = 1.0 • α = 1.5 • (b)
Figure 3 -At left, the spectrum close to the origin of the FSI LSA analyses, for a fixed K dim s = 0.15N m/rad, when varying the angle of attack. At right, the K dim s is varying, such that the equilibrium requirement is met.

reaching the value of x/c = 0.46 for the case at α = 2 • . At this angle, we observe a reattachment of the upper surface separated flow to the rear of the airfoil, giving rise to a recirculation bubble, refereed to as a laminar separation bubble. The aerodynamic loads acting on each x section along the chord are superposed to the velocity field. Each of these infinitesimal forces contributes to the total moment coefficient evaluated around the elastic centre, X EC , whose evolution, with α, is present in figure 2, at right. As depicted by the black arrows in the case α = 1 • , the forces at the rear of the airfoil are pointing downwards, generating an initial negative moment coefficient. As the angle increases, the sense of these forces is reversed, starting from the trailing edge, up to the case α = 2 • , where all contributions have an upward sense, generating a positive moment coefficient. From a static point of view, the initial negative slope of the C mEC -α curve gives us the information that the airfoil has a tendency to pitch-up, as the angle of attack increases. This observation is coherent with the unstable static mode found from the FSI linear stability analysis at α = 0 • . As the angle of attack increases, a second equilibrium can be found, depending on the structural parameters. This second equilibrium position is given by equation 3, after considering a static perturbation in the form θ = θ + θe λt :

       O(1) : K s θ = m z (θ) , O( ) : λ 2 I s θ + λD s θ + K s θ = ∂m z ∂θ θ=θ θ . (8a) (8b)
Noting that α = -θ, equation 8a gives the equilibrium position α, for a given moment coefficient and structure rigidity. The second equation informs us about the presence of an unstable static mode, depending on the slope of the C mEC -α for the considered α. Further, if one considers D s = 0, the stability criterion can be formed as:

           if ∂m z ∂α θ=θ < -K s ⇒ Unstable Static Mode Present , if ∂m z ∂α θ=θ > -K s No Unstable Static Mode Present . (9a) (9b)
The validity of the above stability criterion is first verified without respecting the equilibrium position. For that, an FSI stability analysis is carried out for the different angles with the rigidity parameters I dim s = 0.00135kg m 2 , D dim s = 0 and K dim s = 0.15N m/rad, with the correspondent spectrum present in figure 3a. The results confirm the stability criterion of equation 9, i.e., for the cases α = 0 • and 0.5 • , the slope of the C mEC -α is lower than the rigidity of the airfoil, giving a static unstable eigenmode. On the other hand, for the cases of α = 1 • and 1.5 • , the slope of the C mEC -α higher than the rigidity of the airfoil, resulting in a stable static mode. We note that, although stable, the flutter mode for the latter case has a higher growth rate, being least stable than the lower angle cases.

The same reasoning can be applied by respecting the equilibrium position on equation 8a. For that, and since the values of C mEC and α are fixed by the DNS, the value of K s must be adjusted in order to respect the equilibrium. The results are present in figure 3b, for D s = 0. For the case α = 0.5 • , the mode is expected to be marginally unstable, as the slope of C mEC -α matches exactly the slope rigidity. However, the eigenmode is found to be in the stable side of the spectrum. For the higher cases α = 1.0 • and α = 1.5 • , the static mode becomes stable, as, again, the slope of the C mEC -α is higher than the rigidity of the airfoil. These stable static eigenvalues are not represented in the spectrum as they are mixed with the remaining stable spurious modes. In case at α = 1.5 • , the flutter mode is also found to have moved towards the unstable region, although remaining stable. The change in the frequency of the flutter mode is associated to the fact that the considered rigidity of the structure changes, to match the equilibrium position requirement.

Conclusion

The onset of pitch-oscillations of a NACA0012 airfoil mounted on a torsion-spring was investigated by the means of a global fluid-structure linear stability analysis, around a mean flow field, issued from the time and spanwise averaged DNS, at different incidences, for Re = 50 000. The results at α = 0 • show the presence of a static unstable mode, associated to the divergent departure of the airfoil from its equilibrium position. These results are confirmed by the presence of a negative slope on the C mEC -α, corroborating the tendency of the airfoil to pitch-up. The LSA computations were then carried out for different incidences, taking into account the second equilibrium position attained by the airfoil. As the angle of attack increases, the static mode has become stable, whereas the flutter mode, although stable, has approached the marginal axis. As a future work, the authors suggest to increase the Reynolds number, as the case studied represent the lower boundary of the laminar oscillation flutter instability found by [START_REF] Poirel | Self-sustained aeroelastic oscillations of a NACA0012 airfoil at low-to-moderate Reynolds numbers[END_REF].

  s is the non-dimensional moment of inertia around the elastic centre position, D s the non-dimensional structural damping, K s the non-dimensional structural stiffness and m z the e z -component of the aerodynamic moment, m = m z e z , measured around the elastic centre. The corresponding dimensional structure parameters are noted I dim s

Figure 1 -

 1 Figure1-On the left, the fluid spectrum showing a unstable eigenvalue, with a frequency correspondent to the vortex-shedding frequency. On the centre, a close-up of the same spectrum near the origin. On the right, a fluid-structure spectrum near the origin, highlighting the static and flutter modes, along with the eigenvalue of the structure equation in vacuum, represented by the black cross.

Figure 2 -

 2 Figure 2 -At left, the rear airfoil mean flow topology, illustrated by the streamwise velocity component. The resultant forces from the aerodynamics loads applied to the lower and upper surfaces, that contribute to the moment around X EC , for each section along the chord, are represented by the vertical arrows. At right, the evolution of the moment coefficient of the mean flow with the angle of attack, showing a first range of incidences where the moment is negative.


L Villard

P Borderies

On the use of virtual ground scatterers to localize double and triple bounce scattering mechanisms for bistatic SAR

Keywords: electromagnetic (em) modeling, double and triple bounce scattering mechanisms, image theory, bistatic synthetic aperture Radar (SAR)

In Radar remote sensing, the electromagnetic (em) modeling of scattering contributions from natural or man-made scatterers has often to deal with coupling scattering mechanisms due to the underlying reflecting surface (ground or water). For the monostatic configuration, this contribution is commonly associated with the so-called double bounce scattering mechanism, which is also known to be located in synthetic aperture radar (SAR) images at the vertical projection of the scatterer onto the ground. To explore this phenomenon in the more general case of bistatic SAR but also for triple bounce scattering, a new formalism using virtual scatterers is introduced in this paper. Based upon image theory and on far-field assumption, these virtual points are defined by simple interaction scatterers located on the reflecting surface, and involving the same traveling wave path than the considered scattering mechanism. Depending only on the incidence and scattering angles, elegant formulae in their simplicity are highlighted to construct these fictitious scatterers. The resulting reduction of multiple to simple interaction scattering mechanisms simplifies the finding of their location within bistatic SAR images, generalizing thereby migration effects associated to coupling mechanisms in the bistatic configuration.

Introduction

Electromagnetic (em) scattering from media made of discrete scatterers over an infinite reflecting surface has originated an abundant literature, especially in Radar remote sensing for the modeling of scattering from natural lands (see for instance [START_REF] Ulaby | Microwave remote sensing: active and passive[END_REF]). Considering a discrete medium, the Distorted Born Wave Approximation (DBWA) is commonly applied as long as its volume fraction (volume of scatterers per unit volume within the medium) does not exceed values about 5% (cf. [START_REF] Tsang | Theory of microwave remote sensing[END_REF][START_REF] Tsang | Scattering of electromagnetic waves[END_REF]). In agreement with this one-order approximation, incident or scattered field on or from a scatterer accounts for the direct path together with the averaged contribution of all the scatterers surrounding that path. The most common formulation of this averaged contribution results from the Forward Scattering Theorem as an application of the stationary phase theorem (cf. [START_REF] Ishimaru | Wave propagation and scattering in random media[END_REF][START_REF] Liao | Foldy-lax approximation on multiple scattering by many small scatterers[END_REF]), and leads to the famous 'Foldy' approximation, also referred to as the Effective field approximation. With regard to modeling of natural land surfaces such as agriculture fields or forests, this approximation is broadly used, especially for discrete and coherent (i.e. phase preserving) models such as those detailed in [START_REF] Lin | A Monte Carlo coherent scattering model for forest canopies using fractalgenerated trees[END_REF][START_REF] Marliani | Simulating coherent backscattering from crops during the growing cycle[END_REF][START_REF] Thirion | Capabilities of a forest coherent scattering model applied to radiometry, interferometry, and polarimetry at P-and L-band[END_REF][START_REF] Villard | Forward and inverse modeling of bistatic SAR Observables: Applications in forest remote sensing[END_REF].

To account for the coupling effect with the underlying ground surface -which is modeled as an infinite horizontal plane -the specular point theory (SPT) is commonly used (cf. [START_REF] Beckmann | The scattering of electromagnetic waves from rough surfaces[END_REF]). As described hereinafter, points on which specular reflections onto the ground hold are found using mirroring symmetry with respect to the reflecting surface, these specular points being intrinsic to double and triple bounce scattering mechanisms. At this point, it can be stressed that these restrictions (in the number of interactions accounted for) are not only determined by the growing complexity of these em scattering contributions but also by their isotropy essential for standard SAR compression methods (see e.g. [START_REF] Bamler | Principles of synthetic aperture radar[END_REF]), whether in the common monostatic configuration or in the bistatic one, characterized by a distance between Transmitter and Receiver of the order of their distance to the observed scene (see e.g. [START_REF] Cherniakov | Bistatic radars: emerging technology[END_REF]). For both kind of SAR configurations, the stability of the scattering return for each Radar pulse recorded along the acquisition long or slow-time axis -also called azimuth axis but possibly ambiguous in bistatic geometries -is required in order to get a standard point compression and retrieve the scattering contribution from the target. Hence, the limitations in our study of quasi-isotropic coupling terms (i.e. sufficiently stable to be properly retrieved by standard SAR compression) implying thereby quasi-isotropic scatterers and specular interactions onto the reflecting surface. Higher order and specific coupling terms such as those theoretically detailed in [START_REF] Van Zyl | The Effect of topography on radar scattering from vegetated areas[END_REF] or experimentally analyzed in [START_REF] Auer | Ghost persistent scatterers related to multiple signal reflections[END_REF] will be therefore left aside.

The question that arises in the following lies on the migration of scatterers involved in double and triple bounce scattering in bistatic SAR images. Such question is all the more challenging that severe image distortions may occur in bistatic geometries (cf. [START_REF] Villard | Pol-InSAR simulations in forest bistatic scattering[END_REF]) and that issues related to the so-called ghost scatterers associated with multiple scattering still give rise to many studies in the monostatic case (cf. [START_REF] Auer | Linear signatures in urban sar images -partly misinterpreted?[END_REF]), particularly with respect to their stability as persistent scatterers (e.g. [START_REF] Margarit | Polarimetric characterization and temporal stability analysis of urban target scattering[END_REF]). Given the complexity of multiple scattering, the possibility of finding virtual scatterers involving only one interaction together with the same traveling wave path (and thereby the same phase history) as the original higher order scattering mechanism would thus be very interesting. This is precisely the aim of the following developments, considering on the one hand the image theory framework and on the second hand plane wave fronts rather than spherical ones (hereinafter referred to as far-field approximation).

To introduce the method and its related hypothesis, definitions related to the bistatic configuration are given in the next section with the straightforward case a single reflection off a volume scatterer. Sections 3 and 4 are then dedicated to double and triple bounce interactions, while resulting applications are discussed in the concluding section.

Problem definition considering simple scattering

To introduce the problem and the notations used hereinafter, the straightforward case of a single bistatic reflection -that bounces off a volume scatterer -is considered in this section, as shown in Figure 1. As opposed to ground scatterers, the vertical position of a volume scatterer 'P' is defined by a non-zero height noted 'h'. Ground or volume scatterers are considered as coherent targets, for which scattering contributions seem to come from individual points corresponding to their phase centers. The question at issue lies in finding the position of a virtual ground scatterer which involves the same traveling wave path than the simple bistatic reflection onto the considered volume scatterer.

In the following, a vector from point P to T is described by the symbol -→ PT, distance from P to T by PT, and an unitary vector by û so that | PT|=| -→ PT/PT|=1. Concerning bistatic geometry, this study will be restricted to plane configurations, i.e. for every position Figure 1. Virtual ground scatterer noted E (1) in the case of a single bistatic interaction, considering the em emitted wave for transmitter T which bounces off the volume scatterer P above the ground and then towards receiver R. Assuming that the ellipsoid iso-range ε can be approximated by its tangent, the position of E (1) is given by the intersection between the ground and the tangent driven by the angle ι (1) . of the receiver within the plane defined by ki and ẑ, being, respectively, the unitary incident wave vector and unitary vertical one. The 3D generalization would indeed overwhelm the problem with geometrical features which may be not appropriated to keep the line of this letter. Following that restriction, a bistatic configuration is then completely defined through the plane coordinates of the transmitter and receiver noted 'T' and 'R', respectively. Considering the volume scatterer P (see Figure 1), the incidence and scattering angles θ T and θ R can be emphasized as well as the iso-range loci for the distance |TP+ PR|, which turns out to be the ellipsoid (ε) of foci T and R (cf. [START_REF] Villard | Pol-InSAR simulations in forest bistatic scattering[END_REF][START_REF] Villard | Bistatic Pol-InSAR scenario & evaluation by forest scattering simulations[END_REF]), as the generalization to bistatic configuration of the spherical iso-ranges in the monostatic case. At that stage, it can be noted that the bistatic angle noted and defined by β = (θ Tθ R ) does not fully characterize the bistatic configuration, in particular since a zero bistatic angle may be associated to a non-monostatic geometry, when PT = PR. Besides, whether for airborne or spaceborne configurations, the height noted h of a scatterer above the ground turns out to be sufficiently small in comparison to PT or PR so that the ellipsoid line passing through P can be approximated by its tangent, which intersects with the ground giving the projection into the corresponding pixel of all the points belonging to the same iso-range |TP + PR|. This intersection defines the equivalent ground point for a single interaction onto a volume scatter, noted E (1) . To make its location explicit, one can recall the well-known property of the tangent to an ellipsoid which turns out to be normal to the bisectrix vector β at the tangent point, that is:

d -→ OP dt • PT + P R = d -→ OP dt • β = 0
with t as a real number used to parametrize the vector OP(t). Having in mind this property, it follows that E (1) can be found from the angle ι (1) between the tangent to the ellipsoid and the vertical axis, with:

ι (1) = π 2 -β/2 + θ R = π 2 - θ T + θ R 2 (1)
As a result, this characterization can be considered as a generalization in the bistatic configuration of the well-known layover effect in monostatic SAR images, truly intrinsic to the slant range geometry of Radar acquisitions, unlike optical ones at the nadir of the observed scene which do not suffer from that effect. Indeed, in reducing the 3D scene information through the projection on the ground of all the scatterers belonging to the same iso-range or pixel, the layover effect originates image distortions (e.g. the mixed-up between scattering contributions from several targets) and is thereby essential to account for in monostatic or bistatic SAR imaging.

Virtual ground scatterer for double bounce scattering mechanism

Monostatic case

Even though the monostatic case constitutes a particular bistatic configuration which could have been combined with the next subsection, the developments given hereinafter have not been skipped since it emphasizes the most common case in Radar remote sensing and makes the following geometrical calculations more comprehensive. Within the specular point theory framework (detailed in introduction), the points onto the ground where the specular reflections hold can be judiciously found using the mirroring principle. As shown in Figure 2, considering a given scatterer P and its mirror symmetric noted P ℵ , the specular points are determined by the intersections between the ground surface and the lines (TP ℵ ) and (RP ℵ ) -joining P ℵ with T and R, respectively. Since T and R coincide in the monostatic configuration, the specular points onto the ground noted (G T and G R ) are located at the same position so that the monostatic double bounce is made of two identical paths (with yet opposite time directions), originating thereby coherent effects (superposition of two waves in phase) commonly referred to as the backscattering enhancement (see for instance [START_REF] Saatchi | Coherent effects in microwave backscattering models for forest canopies[END_REF]). In this subsection, T and G T have been arbitrary preferred.

The following two-step demonstration to determine E (2) -the point on which a single interaction involves the same traveling wave path as the double bounce scattering mechanism associated to P -will first consist in showing the existence of E (2) and then its uniqueness. Let us then consider the configuration presented in Figure 2 with E (2) as the vertical projection of P onto the ground. According to the afore-mentioned hypothesis of plane wave fronts (P being sufficiently far for T and R), it can be stated that:

2E 1 E (2) = E 1 E 2 + 2E 2 E (2) + E ′ 2 G T = E ′ 1 P + P ℵ E ′ 2 + E ′ 2 G T = E ′ 1 P + P ℵ G T = E ′ 1 P + PG T ( 2
)
where equality 2E 2 E (2) = P ℵ E ′ 2 results from the Thales' theorem (or also called intercept theorem) applied to the triangle (G T P ℵ P) with (G T P ℵ )//(E 1 E (2) ) and PP ℵ = 2PE (2) , and where equality G T P = G T P ℵ results from construction of P ℵ . From Equation (2), it follows that 2 • TE (2) = TP+ PG T + G T T (3) which proves the existence condition of E (2) . Let's now suppose that equality in (3)i s verified by two points E (2) A and E (2) B . By definition, the simple interaction onto these two points involves the same traveling wave path, so that E (2) A and E (2) B belong to the Figure 2. Monostatic configuration for which transmitter T and receiver R coincide, as well as the specular ground points G T and G R . The virtual ground points for the double and triple bounce scattering mechanisms are, respectively, noted E (2) and E (3) .

same spherical iso-range which intersects with the ground giving a single point, hence the uniqueness condition for point E (2) . Alternatively, it results from the Thales' theorem applied to in the triangles (G T PE ′ 1 ) and (G T P ℵ P) that:

2 • E 1 E 3 = E ′ 1 P and 2 • E 3 E (2) = G T P ℵ = G T P hence: 2 • E 1 E (2) = 2 • E 1 E 3 + 2 • E 3 E (2) = E ′ 1 P + G T P (4)
which also proves the existence condition for E (2) . Finally, a last demonstration can be also put forward since it will be used in the subsequent paragraph in the bistatic case: from the construction of points E 1 and E 4 ,i t results that:

E 1 E 2 = G T E 4 and E 2 E (2) = E 4 P
hence the term by term equality:

G T E 4 + E 4 P + PE ′ 1 = E 1 E (2) + E (2) E 2 + E 2 E 1
which constitutes another necessary condition for the localization of E (2) . Using these developments as an introduction to more advanced cases, the well-known vertical projection on the ground of the double bounce backscattering mechanism has been rewritten. This phenomenon is indeed very common in monostatic SAR imaging: whereas the position of a given scatterer will be shifted according to the layover effect, its double bounce contribution will remain at its vertical projection, whatever the considered scatterer height above the ground.

At that stage, we can already infer that ι (2) -the angle between PE (2) and ẑ -is likely to be proportional to the bistatic angle β, so that β nullifies ι (2) in the monostatic case. (2) corresponds to the position of the floating point E onto the ground so that:

L T → E → R = E=E (2) L T → G T → P → R .
For the sake of clarity, E (2) is not displayed above but can be easily located, moving E on the right so that â = b.

Bistatic case within the plane of incidence

As stated in Section 2, the bistatic configurations considered hereinafter are restricted to the plane of incidence, defined by the incidence wave vector ki and the vertical vector ẑ.The geometry for such configuration is displayed in Figure 3, where an arbitrary point scatterer P above the ground is chosen, as well as transmitter T and receiver R, from now on located at different places. Still considering the specular point theory framework, specular ground points involved for the double bounce scattering mechanism are determined using the mirror symmetric point P ℵ . Two different specular ground points noted G T and G R with respect to T and R can then be found, noticing that both do not coincide anymore as for monostatic cases.

As a result, two distinct equivalent ground points will be searched for, starting from the one with respect to the transmitter which is still noted E (2) in the following (as long as both are not involved). Geometrical features are given in Figure 3, especially for the construction of G T based on P ℵ . Dot lines are dedicated to virtual geometrical points (especially towards E (2) ) whereas plain lines are dedicated to the 'real' traveling wave path: T → G T → P → R, whose length will be hereinafter noted with the operator L • .

Let's consider the angle ι between the vertical axis (along ẑ and the line passing by P and the candidate floating point E, which position on the ground is sought to be E (2) . Points E 1 and E 4 are constructed from the projection of G T and E, respectively, so that whatever the position of E onto the ground, the triangle of base G T E is isosceles and

E 1 E = G T E 4 .
Let us now emphasize the angles â and b, respectively, defined in triangles E 2 EP and EPE 4 (still in Figure 3). These angles can thus be expressed as follows:

â = θ R + ι and b = θ T -ι (5)
With that construction, the equality EE 2 = E 4 P holds if and only if the triangle of base EP is isosceles, i.e. â = b. From that equality, it comes that:

G T P = G T E 4 + E 4 P = E 1 E + EE 2
As a result, with the starting assumption that both transmitter and receiver are far enough from P to consider the equi-phase plane surfaces, a sufficient and necessary condition on the location of E to verify L T → E → R = L T → G T → P → R lies in â = b, hence the following elegant formula in its simplicity for the virtual ground point E (2) associated to the double bounce scattering mechanism with respect to the transmitter, which turns out from ( 5) to be half the bistatic angle:

ι (2) | G T = β 2 (6)
Concerning now the double bounce with respect to the receiver, the reciprocity of the configuration can be invoked: if T and R are inverted, G R is located at the former position of G T so we come back to the traveling wave path associated to the former double bounce with respect to T and then to the former virtual ground point E (2) and its associated angle ι (2) . As expected, reversing the situation by reciprocity does correspond to a switch between G T and G R as well as between T and R in formula [START_REF] Lin | A Monte Carlo coherent scattering model for forest canopies using fractalgenerated trees[END_REF], hence the following generalization:

ι (2) =± β 2 (7)
where the plus or minus are, respectively, for the double bounce with respect to the transmitter or to the receiver. From that result, it is interesting to note that depending on the largest angles between θ T and θ R , the double bounce scattering mechanisms are either projected forward or backward in comparison with the monostatic case where both double bounces are vertically projected on the ground. Besides, this result also means that double bounce contributions (with respect to T or R) of all scatterers located along the line (E (2) P) will be associated with the same traveling distance or time, giving thereby the shapes of Radar isoranges for the scattering mechanisms concerned. Therefore, not only the geometrical phase for bistatic double bounces are different (unlike the monostatic superposition in phase) but these contributions may be not recorded in the same pixel, especially for SAR images at high (spatial) resolution.

Virtual ground scatterer for triple bounce scattering mechanism

The following developments are now dedicated to triple bounce scattering. To come up with the traveling path associated to that scattering mechanism, the coupling interactions between vegetation and ground scatterers are considered for both downwards and upwards paths, respectively from transmitter T towards an arbitrary volume scatterer (still noted P) and from P to receiver R (see Figure 4). Still in agreement with the afore-invoked specular point theory, both specular ground reflections onto points G T and G R are thereby involved, being still constructed using P ℵ , the mirror symmetric of P (as in Section 3). Supported by geometrical details provided in Figure 4, the problem now at issue being in the location of the floating point E onto the ground so that: The virtual ground point E (3) corresponds to the position of the floating point E onto the ground so that:

L T → E → R = E=E (3) L T → G T → P → G R → R (8)
L T → E → R = E=E (3) L T → G T → P → G R → R .
In agreement with the previous hypothesis, equi-phase plane surfaces are still considered so that:

L T →E → R -L T → G T → P → G R → R = E 1 E + EE 2 -G T P + PG R = E 1 E 3 + E 3 E + E ′ 4 G R -G T P + PG R = E 1 E 3 + E 3 E + E ′ 4 G R -G T P ℵ + P ℵ G R = E 3 E + E ′ 4 G R -P ℵ G R = E 3 E -P ℵ E ′ 4 (9)
From the above equation, it can be deduced from the last line that equality [START_REF] Thirion | Capabilities of a forest coherent scattering model applied to radiometry, interferometry, and polarimetry at P-and L-band[END_REF] holds if

E 3 E = P ℵ E ′ 4 , i.e. if the line P ℵ E bisects the corner E ′ 4 P ℵ E 4 , meaning that: a = b so that: θ R + ι = π -θ T -ι (10)
providing the sought formula:

ι (3) = θ T + θ R 2 - π 2 =-ι (1) (11)
Existence condition for E (3) is thus proven by the above expression. In the event that two points E (3)

A and E

(3) B exist, these points onto which simple interactions occur would belong to the same ellipsoid iso-range, which intersects with the ground giving a single point, proving thereby the uniqueness condition for point E (3) .

As far as SAR applications are concerned, the magnitude of triple bounce scattering mechanism is more likely to be significant over built-up structures such as ships, bridges, or buildings (see for instance [START_REF] Auer | Linear signatures in urban sar images -partly misinterpreted?[END_REF][START_REF] Soergel | Feature extraction and visualization of bridges over water from high-resolution insar data and one orthophoto[END_REF]) as well as flooded regions, where specular ground reflections are favored by flat and smooth surfaces.According to the previous result Equation [START_REF] Bamler | Principles of synthetic aperture radar[END_REF], it is rather remarkable that whatever the scatterer above the ground, its triple bounce contribution will appear backwards, with exactly the opposite distance which shifts forward the volume contribution, generalizing thereby the layover effect to triple bounce scattering mechanism in bistatic configuration.

Discussion on resulting applications & conclusion

In this paper, an elegant theoretical method to reduce double and triple bounce scattering mechanisms to a simple interaction onto a virtual ground scatterer is developed, for any bistatic configuration where the receiver belongs to the plane of incidence. This method is based on virtual scatterers located on the ground surface and for which the simple interaction involves the same traveling wave path as its associated multiple scattering mechanism. With respect to the original scatterer, the location on the ground of this virtual point turns out to be driven through remarkable formulae in their simplicity, involving simply half the signed bistatic angle in the case of double bounce scattering mechanismwith opposite sign between double bounce with respect to transmitter or to receiver -and the opposite angle associated with the simple layover in the case of triple bounce scattering mechanism.

Beyond novelty of the theoretical formulation in SAR configurations for which we have little experience, these results offer a simple but efficient way to characterize the migration phenomenon in bistatic images and the so-called ghost scatterers, mostly due to multiple scattering. However, as for current state-of-the art monostatic SAR imaging methods, non-supervised corrections to avoid these ghost scatterers are not applicable without a-priori information about the observed scene. In addition, these results also demonstrate that the voxel shapes (i.e. the 3D volumes which encompass all the scatterers whose specific contribution falls into a given pixel) is intrinsic to each scattering mechanism, which is essential for analytical formulation in em modeling and especially to derive the interferometric phase sensitivity to scatterer height.

Being less and less limited by technological barriers (especially related to synchronisation), one can currently attest of a new resurgence of bistatic Radar, not only characterized by airborne concepts but also by hybrid configurations (like GNSS-Reflectometry) and the willing of putting bistatic configurations into space (e.g. with the SAOCOM-CS ESA candidate mission [START_REF] Gebert | SAOCOM-CS -a passive companion to saocom for single-pass L-band sar interferometry[END_REF]). Fostered by these promising concepts, next steps will be dedicated to generalize our results to the full bistatic hemispherical space (combining geometrical features of both scattering and incident planes) and to specific SAR applications, in the case of non-stationary configurations and polarimetric and interferometric bistatic SAR acquisitions.

Figure 3 .

 3 Figure 3. Traveling wave path in the bistatic configuration for the double bounce scattering mechanism with respect to the transmitter T. The virtual ground point E(2) corresponds to the position of the floating point E onto the ground so that: L T → E → R =

Figure 4 .

 4 Figure 4. Traveling wave path in the bistatic configuration for the triple bounce scattering mechanism.The virtual ground point E(3) corresponds to the position of the floating point E onto the ground so that: L T → E → R =

Acknowledgements

The authors are truly grateful to Thuy Le Toan, especially for its interest and support which ease the finalization of developments and writing. This study has been also made possible thanks to the financial support of CNES (Centre National d'Etudes Spatiales) through the TOSCA research program.


Edison Gerena
email: edison.gerena@isir.upmc.fr

Stéphane Régnier

Sinan Haliyo
email: sinan.haliyo@sorbonne-universite.fr

Interactive Laser-actuated micro-robots for Experimental Biology

Recent breakthroughs in biotechnology are rising the demand for single-cell manipulation techniques such as cell isolation and 3D orientation. Significant challenges remain for these applications, mainly due to the physics involved and the size constraints imposed by the environment. The precision required at those dimensions have a cost in term of workspace, grasping strategies and control schemes. Furthermore, there is an increasing demand to manipulate objects in confined environments like micro-fluidic devices, rendering external actuators unusable. In this presentation, an interactive optical robotic system is introduced. It implements unthetered actuation and wireless sensing for microobjects in a confined liquid medium.

Introduction

Cell Biophysics is an emerging interdisciplinary science, in which physics and chemistry are employed to understand biological mechanisms, to assay cell response to events, characterize activation status in different settings or coming from different healthy or patient individuals. As such, adding to established spectroscopic and micro-scopic techniques, a new range of mechanical techniques and devices has emerged in recent years to investi-gate, measure and map mechanical forces in cells and decipher their role in cell function and fate. In this context, the ability to manipulate and characterize individual biological objects is become a major scientific challenge. The need to perform tasks such as single-cell deformation, stimulation, rotation, or transportation, has called for new emerging techniques that allow working at the unicellular scale.

Among these devices, robotic micromanipulators have become an ubiquitous tool in many experiments as they offer the only means to access the samples for manipulation, physical characterization, extraction… Using needles or glass pipettes as end-effectors, robotic systems give the opportunity to directly process a sample such as an isolated cell, a bacterium, making possible a wide range of mechanical, chemical or electrical stimuli and measurements [START_REF] Chaillet | Microrobotics for Micromanipulation[END_REF]. However, these tools are complex to integrate in an experimental set-up and have dimensions many times larger than samples. It is thus quite complicated to design such an experiment as well as mastering the proper use of these devices. Additionally, using electro-mechanical sensors directly on the samples is quite difficult as those are prone to be affected by the liquid environment, provide low sensitivity and precision, and low signal-tonoise ratio. In addition, interfacing and communicating with the physical sensors through the sample chamber is highly challenging.

Optical techniques may greatly overcome these difficulties by providing wirelessly both the actuation and sensing. Optical tweezers [START_REF] Ashkin | Acceleration and trapping of particles by radiation pressure[END_REF] avoid many of the limitations of competing techniques as they offer much greater versatility. Optical trapping relies on an immaterial electromagnetic field produced by a highly focused laser beam. A highly focused laser produces a localized three-dimensional electromagnetic field that stably traps microscopic objects. This object is then actuated by deflecting the focal point of the laser. The optical forces can be modeled, and three-dimensional trap stiffness can be estimated analytically or numerically [START_REF] Neuman | Optical trapping[END_REF]. Motion tracking allows the force to be inferred. A trapped microsphere can serve as a force probe to examine the physical characteristics of a target object, such as the stiffness of cells, thus enabling a broad variety of potential applications, especially in biology or medical science. Optical Tweezers have successfully been applied in a large range of bio-manipulation experiments [START_REF] Zhang | Optical tweezers for single cells[END_REF], such as the trapping of red blood cells and the immobilization of bacterial cells for nanoscopy.

Optical Microrobots is an evolution of optically trapped simple probes mentioned above. These are microscale floating 3D structures elaborated with optical 'handles' and actuated by optical tweezers. Using the laser tweezers as the actuation principle, and real-time motion tracking for sensing, these can be controlled in closed-loop or be teleoperated, in complete 6 degrees of freedom.

Interactive optical robotic system

In order to take the full advantage of these optical robots, our team has produced a robotic optical trapping system by carefully take in to account the sensor feedback, actuation and control issues (Fig. 1.A). It include a realtime 3D force measurements using event-based vision sensors. These bio-inspired sensors are frame-free and eliminate data redundancy by design and allow for 2D tracking at the rate in the order of tens of kilo-Hertz [START_REF] Ni | Asynchronous event-based visual shape tracking for stable haptic feedback in microrobotics[END_REF]. The integration of such a sensor in laser trapping set-up provide high-bandwidth 3D motion tracking [START_REF] Yin | High-bandwidth 3D force feedback optical tweezers for interactive bio-manipulation[END_REF]. External forces applied on the trapped object are inferred with piconewton resolution in nanonewton range in practically real-time, up-to 300 Hz sampling rate. Haptic exploration tasks in biological (Fig. 1.B) and not-biological objects [START_REF] Gerena | Improving Optical Micromanipulation with Force-Feedback Bilateral Coupling[END_REF] has been performed.

An active solution for 3D high speed actuation of optical tweezers with MOEMS [START_REF] Gerena | High-bandwidth 3-D multitrap actuation technique for 6-DoF real-time control of optical robots[END_REF] has also be proposed. It allows to generate more than 15 time-shared traps in 3D with low latency and high bandwidth. The 3D motion of the focal spot is obtained by the synchronization of the orientation of a galvanometer mirror and the focusing or defocusing of a deformable mirror. This actuation technique is useful in applications where the 3D orientation of microscopic objects is needed, such as cell surgery or 3D tomographic imaging of living samples. The individual 3D manipulation of multiple objects, and the 6DoF manipulation of optical robots is depicted in Fig 1 .C. A teleoperation interface for interactive bio-manipulation is also implemented [START_REF] Gerena | Tele-Robotic Platform for Dexterous Optical Single-Cell Manipulation[END_REF]. The system provides a straightforward human/machine interaction through a tele-robotic solution allowing dexterous manipulation of synthetic and biological objects in an efficient and intuitive way. Teleoperation control is implemented with an Omega 7, a haptic device with 7-DoF (including an active two-finger grip) and force feedback, and is ensured by a hard real-time operating system. Traps can be grouped and controlled in a variety of ways for specific purposes.

Optical robots have dimensions in the same range of samples and can be simply introduced in the medium. Two diferents examples are show in Fig. 1.D. Compared to simple spherical probes, these can be equipped with custom designed end-effectors, for example sharp probes, mechanical sensors or injectors, shovels... to carry out different stimuli and measurements directly on the bio-samples. These Optical robots can be used collectively to induce mechanical stimuli, squeeze, or bring samples in contact and directly measure their responses at high temporal and force resolution, with eventually through haptic devices for "feeling" the force and dynamics of the cell surface. Looking further forward, they can include internal mobilities, to do tools like grippers, or active components like a micromotor or a pump. This represents a huge potential for experimental biology, as the user will virtually be able to reach inside the Petri Dish and manipulate single specimens, while collecting information on interaction forces.


Bounds on Herglotz functions and physical limits to broadband passive cloaking in the quasitatic regime Maxence Cassier (speaker) (joint work with Graeme W. Milton) Introduction: Cloaks are specific structures placed near or around an object that render the electromagnetic response of cloak plus object equal or almost equal to that of free space. Ideally passive cloaks should work for waves in a broadband frequency range, giving rise to the challenging question: Is it possible to perform broadband passive cloaking over a finite frequency band? In the context of the quasistatic approximation of Maxwell's equations we prove that it is impossible and give quantitative limitations to cloaking over a finite frequency range. Our results, published in [START_REF] Cassier | Bounds on Herglotz functions and fundamental limits of broadband passive quasistatic cloaking[END_REF], hold for a cloak or object of any geometrical shape and do not depend on the cloaking methods: transformation optics, anomalous resonance, complementary media.

1) The passive cloaking problem

Let O be a bounded simply-connected dielectric inclusion with Lipschitz boundary that one wants to cloak. O is characterized by its permittivity ε(x, ω) = ε I, where ε > ε 0 is constant on the frequency range of interest [ω -, ω + ] and strictly larger than the permittivity of the vacuum ε 0 . The passive cloak is made of an anisotropic material of any shape characterized by its dielectric tensor ε(x, ω) which depends both on the spatial variable x and the frequency ω. The whole device, the inclusion and the cloak, occupies an open bounded set Ω ⊂ B(0, R 0 ) of characteristic size R 0 and the remainder of space R 3 \ Ω is vacuum of permittivity ε(x, ω) = ε 0 I. The observer is assumed to be at a distance R R 0 . We send a plane wave towards the device and assume that its wavelength is considerably larger than R in the frequency range of interest ω ∈ [ω -, ω + ] ⊂ R +, * so that we can use the quasistatic approximation in this frequency band. In this approximation, the curl-free electrical field E(x, ω) is given in terms of the gradient of some potential V (x, ω), i.e. E(x, ω) = -∇V (x, ω), the incident plane wave in the vicinity of a closed ball B(0, R) corresponds to a uniform field E 0 ∈ C 3 so that the potential ∇V (x, ω) satisfies the following elliptic equation

∇ • ε(x, ω)∇V (x, ω) = 0 on R 3 , (1)
and admits the asymptotic expansion as |x| → ∞:

V (x, ω) = -E 0 • x + p(ω) • x 4πε 0 |x| 3 + O 1 |x| 3 , with p(ω) = α(ω)E 0 . (2)
Thus, the main contribution of the scattered far field is a dipolar term p(ω) ∈ C 3 which depends linearly on E 0 via the polarizability tensor α(ω) ∈ M 3 (C). Hence, to cloak the device Ω at a sufficient large distance R to any incident field E 0 ∈ C 3 at a frequency ω ∈ [ω -, ω + ], one needs that α(ω) vanishes at ω. The electric induction D is given within the cloak Ω \ O by the constitutive law (CL): D = ε 0 E+ε 0 χ E t E, where t stands for the time convolution product 1 between the real-valued susceptibility tensor χ E (x, t) and the electrical field E. To define this convolution, one assumes for simplicity that

χ E ∈ L 1 (R t , L ∞ (Ω \ O) 9 ) and that E, ∂ t E ∈ L 2 (R t , L 2 (Ω \ O) 3 ). The cloak is a passive material since it is causal: χ E is supported in (Ω \ O) × R +
t and passive, i.e. for any real fields (E, D) satisfying the (CL) (and the regularity assumptions for E):

t -∞ Ω\O ∂ t D(x, t) • E(x, t)dxdt ≥ 0, ∀t ∈ R. Let C + := {ω ∈ C | Im ω > 0} and cl C + := C + ∪ R. For any causal f ∈ L 1 (R t ),
one defines the Fourier-Laplace transform as f (ω) := R + f (t)e iωt dt, ∀ω ∈ cl C + so that it coincides with the Fourier transform for real frequency. The (CL) in the frequency domain becomes: D(x, ω) = ε(x, ω) Ê(x, ω) with ε(x, ω) = ε 0 (1 + χE (x, ω)), ∀ω ∈ R. Thus, one shows that the passivity of the cloak is equivalent in the frequency domain to ( H1 ): for a.e. x ∈ Ω \ O, ε(x, •) is analytic on C + and continuous on cl C + , ( H2 ): for a.e.

x ∈ Ω \ O, ∀ω ∈ cl C + , ε(x, -ω) = ε(x, ω), ( H3 ): for a.e. x ∈ Ω \ O, ∀ω ∈ R + , Im ε(x, ω) ≥ 0 (passivity), ( H4 ): for a.e. x ∈ Ω \ O, ε(x, ω) → ε 0 I as |ω| → ∞ in cl C + .

2) Bounds on Herglotz functions and passive systems

Herglotz functions are analytic functions of the upper-half plane with nonnegative imaginary part. In [START_REF] Cassier | Bounds on Herglotz functions and fundamental limits of broadband passive quasistatic cloaking[END_REF], we derive bounds on Herglotz functions which apply to a wide class of linear passive systems and generalize those provided in [START_REF] Bernland | Sum rules and constraints on passive systems[END_REF]. To this aim, we consider a passive linear system characterized by a function f : cl C + → C in the frequency domain which satisfies the assumptions

(H 1 ) f is analytic on C + , continuous on cl C + , (H 2 ) f (-ω) = f (ω), ∀ω ∈ cl C + , (H 3 ) Im f (ω) ≥ 0, ∀ω ∈ R + , (H 4 ) f (ω) → f ∞ > 0, when |ω| → ∞ in cl C + ,
i.e. ( H1 -H4 ) but for a scalar function. We define the square root by √ ω = |ω| 1 2 e i arg ω/2 if arg ω ∈ (0, 2π) and by

√ x = |x| 1 2 for x ∈ R + . In [2], we show that v defined by v(ω) := ωf ( √ ω), ∀ω ∈ C is a Herglotz function that is analytic in C \ R + , negative on R -, * and satisfies by (H 2 ): v(ω) = v(ω), ∀ω ∈ C + ∪ R -, * and by (H 4 ): v(ω) = f ∞ ω + o(ω) when |ω| → ∞ in C + .
Then, we introduce the Herglotz functions h m and v m defined by:

h m (ω) = R dm(ξ) ξ -ω and v m (ω) = h m (v(ω)), ∀ω ∈ C + ,
where m ∈ M with M is the set of probability measures on R. Using a sum rule derived in [START_REF] Bernland | Sum rules and constraints on passive systems[END_REF], we show (see [START_REF] Cassier | Bounds on Herglotz functions and fundamental limits of broadband passive quasistatic cloaking[END_REF]) the following theorem.

Theorem 1. Let [x -, x + ] be a compact interval of R + * , then one has:

lim y→0 + 1 π x+ x- Im v m (x + iy) dx ≤ 1 f ∞ ∀m ∈ M, (3)
and Dirac measures (δ ξ ) ξ∈R optimize the inequality (3) since

sup m∈M 1 π lim y→0 + x+ x- Im v m (x + iy) dx = sup ξ∈R 1 π lim y→0 + x+ x- Im v δ ξ (x + iy) dx.
We prove this result in [START_REF] Cassier | Bounds on Herglotz functions and fundamental limits of broadband passive quasistatic cloaking[END_REF] for compactly supported measures m ∈ M but it can be shown exactly in the same way for any m ∈ M. If Im f (ω) = 0 for ω ∈ [ω -, ω + ], we show by using Dirac measures (δ ξ ) ξ∈R in inequality (3) that

ω 2 0 (f (ω 0 ) -f ∞ ) ≤ ω 2 (f (ω) -f ∞ ), ∀ω, ω 0 ∈ [ω -, ω + ] such that ω 0 ≤ ω. (4)
Without such assumption on Im f , by using the uniform probability measure on

[-∆, ∆] with ∆ = max [ω 2 -,ω 2 + ] |v(x)| in the bound (3) it follows that: 1 4 (ω 2 + -ω 2 -)f ∞ ≤ max ω∈[ω-,ω+] |ω 2 f (ω)|. (5)
.

3) Fundamentals limits to broadband cloaking

We apply now the above bounds to cloaking. We give in [START_REF] Cassier | Bounds on Herglotz functions and fundamental limits of broadband passive quasistatic cloaking[END_REF] a functional framework to equations ( 1) and ( 2) (which are physically relevant in [ω -, ω + ] where the quasistatic approximation is valid but holds for any ω ∈ cl C + by using the analytic extension of the permittivity in the inclusion and in the vacuum).

Then, we show (with a coercivity assumption, see [START_REF] Cassier | Bounds on Herglotz functions and fundamental limits of broadband passive quasistatic cloaking[END_REF]) that for a passive cloak satisfying ( H1 -H4 ) and a reciprocity principle, the function f E0 given by

f E0 (ω) := α(ω)E 0 • E 0 = Ω (ε(x, ω) -ε 0 I) E(x, ω) • E 0 dx, ∀ω ∈ cl C + is well-defined for E 0 ∈ C 3 and satisfies ( H1 -H4 ) with f E0,∞ := α(∞)E 0 • E 0 where α(∞) := lim |ω|→+∞ α(ω)
is the positive definite polarizability tensor of the inclusion O depending only on its geometry and contrast in permittivity. If the cloak is a lossless (i.e if Im ε(x, ω) = 0 inside the cloak on [ω -, ω + ]), one shows that Im f E0 (ω) = 0 on [ω -, ω + ]. Thus, using (4) on functions f E0 gives

ω 2 0 (α(ω 0 ) -α(∞)) ≤ ω 2 (α(ω) -α(∞))
, ∀ω, ω 0 ∈ [ω -, ω + ] such that ω 0 ≤ ω, which turns to be an optimal bound (see [START_REF] Cassier | Bounds on Herglotz functions and fundamental limits of broadband passive quasistatic cloaking[END_REF]). Now assume that one can cloak at a frequency ω 0 . Thus α(ω 0 ) = 0 and it yields to α(ω) ≤ -α(∞) ω 2 0 -ω 2 ω 2 , ω ∈ [ω -, ω 0 ] and α(∞)

ω 2 -ω 2 0 ω 2 ≤ α(ω), ω ∈ [ω 0 , ω + ]
which obviously forces α(ω) to be non-zero away from the frequency ω 0 on [ω -, ω + ] and makes cloaking impossible on [ω -, ω + ]. If the cloak is not lossless, one applies the bound (5) on the functions f E0 to get:

1 4 (ω 2 + -ω 2 -)α(∞)E 0 • E 0 ≤ max ω∈[ω-,ω+] |ω 2 α(ω)E 0 • E 0 |, ∀E 0 ∈ C 3 .
This positive lower bound also gives a limitation to cloaking on [ω -, ω + ].


Thomas Hermans

Nicolas Compaire

Robin Thibaut

Nolwenn Lesparre

Bayesian Evidential Learning: an alternative to hydrogeophysical coupled inversion

Deterministic geophysical inversion suffers from a lack of realism because of the regularization, while stochastic inversion allowing for uncertainty quantification is computationally expensive. In this contribution, we propose to use Bayesian Evidential Learning as an alternative to hydrogeophysical coupled inversion. We demonstrate the ability of the approach to successfully predict a hydrogeological target from time-lapse ERT data in the context of a heat injection and storage experiment.

Introduction

The current paradigm in geophysical imaging relies on inverse modelling yielding the "best" solution reproducing the data within its estimated noise level. Since geophysical inverse problems are ill-posed, and solutions are non-unique, the most commonly used approach is to solve a regularized version of the problem, where the data misfit is minimized under some model constraints (Tikhonov and Arsenin, 1977). The image produced is generally the result of a nonlinear inversion or process and represents only one possibility amongst an infinite number of solutions (Linde et al., 2015). Recently, the nature and extent of the non-unicity of the solution, related to the intrinsic uncertainty of the subsurface and the noisy nature of the data, has been acknowledged. The need to account for uncertainty at each step of the imaging process has been recognized (Linde et al., 2015). Therefore, the development of stochastic methodologies for the inversion of geophysical data has become mandatory (e.g., Wathelet, 2008;Trainor-Guitton and Hoversten, 2011;Lochbühler et al., 2014). The currently most widespread approach to stochastically solving the geophysical inverse problem is to use Markov chain Monte Carlo (McMC) methods (Sambridge and Mosegaard, 2002;Vrugt et al., 2013). McMC methods are based on algorithms sampling the posterior distribution of the model parameter space. However, these methods are very computationally demanding, because they typically require hundreds of thousands of forward model runs based on a chain of models dependent on each other. Each step in the process requires to solve the forward problem to estimate the likelihood of the current model (i.e., its data misfit), and hundreds of thousands of iterations might be needed to converge, especially for complex subsurface models (e.g., Irving and Singha, 2010).

However, obtaining a single image is often not the end-goal of a geophysical study. Instead, a geologist, hydrogeologist or reservoir-engineer often wants to derive some interesting (hydro)geological features from it. Post-processing techniques are then used to extract features from this image. Interpretation often relies on uncertain empirical petrophysical relationship (Rubin and Hubbard, 2005;Paasche et al. 2006) and suffers from the spatially-dependent resolution of the image (Day-Lewis et al., 2005). Taking this into consideration when interpreting a geophysical image generally requires the re-introduction of the concept of uncertainty or stochasticity (e.g., Moysey et al., 2005;Singha and Moysey, 2006;Hermans and Irving, 2017) The spatially-dependent resolution of geophysical images and the difficulty to extract quantitative features out of it gave birth to coupled inversion approaches, where geophysical data are integrated to hydrogeological (or other dynamical) inverse problems through a petrophysical relationship (Irving and Singha, 2010;Christensen et al., 2016). Such an approach allows to produce subsurface models consistent with both geophysical and hydrogeological data and to integrate prior geological information. However, parameterizing and calibrating such a model remains very complex and coupled hydrogeophysical inversion often require simplifying assumptions (Christensen et al., 2016) or require computationally expensive stochastic approaches (Irving and Singha, 2010;Hermans et al., 2015).

In this contribution, we propose to circumvent coupled hydrogeophysical inversion through Bayesian Evidential Learning (BEL;Hermans et al., 2016Hermans et al., , 2018)). The rationale behind this approach is that model inversion is rarely the final objective of modelling but rather a mean to generate Table 1. Parameters of the prior model some predictions related to the study site. In BEL, a statistical relationship is sought between data and prediction through physically-based forward modelling, using a limited number (typically a few hundreds) of realizations sampled from the prior model space. This allows to maintain the geological realism of the model and prediction, while replacing the complex inversion in the model space by a simple regression in the joint data-prediction space. This allows to estimate the posterior distribution of the prediction at limited costs. In this contribution, we illustrate the method to predict the thermal affected zone during a heat injection and storage experiment carried out in an alluvial aquifer (Lesparre et al., 2019).

Method

We refer to Hermans et al. (2016Hermans et al. ( , 2018) ) for a detailed description of BEL applied to geophysical data. Here, we only shortly describe the different steps (Figure 1):

1) The prior model is defined, i.e., a prior uncertainty range is given for all model parameters. In our cases, the hydraulic conductivity field is simulated by sequential Gaussian simulation, with uncertain mean, standard deviation and variogram range, anisotropy and orientation. In addition, the amplitude of the natural gradient and the effective porosity are also uncertain (Table 1). We sample 250 realizations from the prior. 2) The experiment is simulated for each sample. It consists in the injection of heated water (Δ𝑇 = 28.6 𝐾) at a rate of 3 m³/h during 6 hours at 5 m depth in a 10 m thick aquifer (Figure 2), followed by a storage period (91 h), a pumping period (15.5 h at 3 m³/h) and a final resting period (Lesparre et al., 2019). Note that if data and prediction correspond to different experiments, both should be simulated.

3) The temperature distribution is extracted from the heat transfer simulation at every time step to define the prediction (thermal affected zone) and the data (time-lapse ERT) collected using 6 parallel profiles of 21 electrodes spaced by 2 m (Figure 2). We use a linear relationship for the dependence of the electrical conductivity to the temperature (Hermans et al., 2014). 4) The original data set contains 1948 quadrupoles, and the prediction is the temperature in a volume composed of 3808 cells, both for 106 time-steps. We first reduce the dimensions of both variables using principal component analysis (PCA). 99.9% of the data and 91.3 % of the prediction are explained by 30 dimensions. 5) We apply canonical correlation analysis (CCA) (Figure 3). This process yields a set of independent linearized relationship between the reduced dimensions of the data (d) and the prediction (h). 6) We sample the joint data-prediction distribution in the reduced dimension space by applying the same dimension reduction to the observed data. Noise is propagated using a Monte Carlo approach to estimate the covariance matrix in the reduced space (Hermans et al., 2016). In contrast to Hermans et al. (2016Hermans et al. ( , 2018)), we do not use a linear regression because the distributions are not Gaussian. Instead, we use Kernel Density Estimation (Hermans et al., 2019, Michel et al., 2020). This provides the posterior distribution of the prediction in the reduced space. 7) The posterior distribution is back-transformed in the original physical space.

Results

The validity of the approach is demonstrated using a synthetic case. One of the realizations of the prior is taken as the true model and the corresponding data set as the observed field data set.

Figure 4 shows a comparison between the true prediction, the median posterior prediction and a deterministic inversion using the standard smoothness constraint. We can clearly observe that samples from the posterior distribution are In particular, the hot plume is vertically limited to the thickness of the screen interval. In contrast, the smoothness constraint inversion underestimates the amplitude of the maximum temperature and spreads out the thermal affected zone over a much larger volume, illustrating the typical drawback of deterministic inversions.

BEL also allows to estimate the uncertainty range, as it is easy to sample any number of samples from the posterior distribution. Figure 5 shows the temporal evolution of the temperature in the vicinity of the injection borehole for the prior realizations, the true prediction and the 5%-95% confidence interval deduced from BEL posterior distribution. We see that the median posterior sample is close to the true model and that the uncertainty range is maximum et the end of the injection phase and in the storage phase, with a value about 1°C. Note that 1°C correspond to a change of electrical resistivity of about 2%. The uncertainty range is much narrower during the pumping phase and for the final resting period.

We estimated the average estimation error of BEL through time (Figure 5, right). The percentage on the figure indicates the number of time steps for which the true prediction is not in the 5-95% interval. Blank cells are cells that are always within this range. We see that deviations occur mostly for grid cells located outside of the center of the plume (Figure 3), i.e., for grid cells that are not experiencing large changes of temperature. In other words, those cells correspond to zones with a very small variation amplitude, which is beyond the resolution of ERT (< 0.5-1 °C). The few cells located in the middle of the model close to the injection borehole are due to the slight overestimation tendency after the pumping phase (Figure 5, left). In comparison, the smoothness constraint inversion lies outside of the confidence interval for every time step.

Conclusion

In this contribution, we propose to use Bayesian Evidential Learning as an alternative to hydrogeophysical inversion. Instead of inverting the geophysical data set to get the posterior distribution of the model parameters, which would be extremely difficult and computationally expensive, we circumvent the inversion by learning a direct relationship between the geophysical data (time-lapse ERT) and the hydrogeological prediction (temperature distribution). It is then straightforward to generate the posterior distribution of the prediction for the observed data set.

This approach is particularly efficient. We only need 250 samples out of the prior distribution to learn the statistical relationship. Therefore, the posterior distribution can be calculated by running only 250 times the forward heat transport model. This is several orders of magnitude less that what would be necessary for an McMC approach. This is possible because the prediction is much less complex than the model parameter distribution. The latter is a heterogeneous hydraulic conductivity distribution with a complex spatial correlation (see Table 1), while the temperature distribution is relatively smooth and of limited volume because of the conduction processes taking place in the subsurface (heat exchange between the solid matrix and the fluid).

The framework is easily adaptable, so not only can it be applied to hydrogeological data (Hermans et al., 2019), but also could be extended to combine geophysical and hydrogeological data. It can also be applied to prediction corresponding to another experiment or even another location in the aquifer (Hermans et al., 2018). Therefore, BEL is a suitable framework to replace hydrogeophysical inversion when the distribution of the model parameters is not an objective of the study.

Figure 1 :

 1 Figure 1: Illustration of the different steps of the BEL framework

Figure 2 :

 2 Figure 2: Experimental set-up

Figure 3 :

 3 Figure 3: Results of canonical correlation analysis (16 first components)

Figure 4 :

 4 Figure 4: This caption is placed inside the frame

Figure 5 :

 5 Figure 5: Estimation of the average temperature around the injection well throughout the experiment (left) and distribution of the average error of estimation (right).