Philippe Laurenc ¸ot

Christoph Walker
email: walker@ifam.uni-hannover.de

VARIATIONAL SOLUTIONS TO AN EVOLUTION MODEL FOR MEMS WITH HETEROGENEOUS DIELECTRIC PROPERTIES

Keywords: October 8, 2019. 2010 Mathematics Subject Classification. 35K86, 74H20, 35Q74, 35M86, 35K25 MEMS, gradient flow, transmission problem, obstacle problem, fourth-order equation. Partially

published or not. The documents may come

Variational solutions to an evolution model for MEMS with heterogeneous dielectric properties

Philippe Laurençot, Christoph Walker

INTRODUCTION

The existence of variational solutions is shown for an evolution problem describing the spacetime dynamics of a microelectromechanical system (MEMS) with heterogeneous dielectric properties. Specifically, a MEMS device such a switch is made of a thin rigid conducting plate above which a thin conducting elastic plate is suspended and clamped at its boundary. The shape of the undeformed elastic plate is identical to that of of the rigid one. Holding the two plates at different electrostatic potentials generates a deformation of the top elastic plate to compensate the induced electrostatic force. It is by now well-known that a sufficiently large potential difference can lead to a pull-in instability or touchdown, a situation which corresponds to the top plate coming into contact with the bottom one and results in a short circuit due to the potential difference [START_REF] Ambati | Some studies on the deformation of the membrane in an RF MEMS switch[END_REF][START_REF] Bernstein | Modeling and analysis of hysteresis phenomena in electrostatic zipper actuators[END_REF][START_REF] Ch | Some singular equations modeling MEMS[END_REF][START_REF] Pelesko | Mathematical modeling of electrostatic MEMS with tailored dielectric properties[END_REF][START_REF] Pelesko | Modeling MEMS and NEMS[END_REF]. Clearly, such a phenomenon may alter the properties or the operating conditions of the MEMS device. However, it can be prevented, for instance, by covering the ground plate with an insulating layer of positive thickness [START_REF] Ambati | Some studies on the deformation of the membrane in an RF MEMS switch[END_REF][START_REF] Bernstein | Modeling and analysis of hysteresis phenomena in electrostatic zipper actuators[END_REF][START_REF] Lindsay | Regularized model of post-touchdown configurations in electrostatic MEMS: Equilibrium analysis[END_REF][START_REF]Regularized model of post-touchdown configurations in electrostatic MEMS: Interface dynamics[END_REF]. We consider this situation herein and thus assume that the bottom plate is covered by a non-deformable layer of positive thickness, possibly having heterogeneous dielectric properties characterized by a permittivity σ 1 > 0, which differs from the constant permittivity σ 2 > 0 of the surrounding medium. Touchdown may still occur, in the sense that the top plate may come into contact with the upper side of the insulating layer. However, such a situation does not generate a singularity, as the top plate cannot penetrate the layer. Assuming further that the physical state of the MEMS device is fully described by the deformation u of the top plate and the electrostatic potential ψ u between the two plates, the dynamics of the MEMS is then governed by the competition between mechanical and electrostatic forces, and is given by a time relaxation towards critical points of the total energy, the latter including mechanical, contact, and electrostatic contributions.

To convert this rough description of the model into mathematical equations, we assume that there is no variation in one of the two horizontal directions and consider a two-dimensional MEMS device in which the rigid ground plate and the undeformed elastic plate have the same one-dimensional shape D := (-L, L), L > 0, the former being located at height z = -Hd, v

Ω 1 Ω 1 Ω 2 (v) D Σ z -H -d -H 0 2L G(v) FIGURE 1.
Geometry of Ω(u) for a state u = v with empty coincidence set (green).

H > 0, d > 0, while the latter is clamped at its boundary (±L, 0). The bottom plate is then D × {-H -d} and it is covered by an insulating layer

Ω 1 := D × (-H -d, -H)
of positive thickness d. The dielectric properties of Ω 1 are characterized by the permittivity σ 1 > 0 which may vary with the horizontal coordinate x ∈ D and the height z ∈ (-Hd, -H). As for the elastic top plate, we assume that the relevant physical framework is restricted to small deformations and, denoting the vertical deformation of the top plate by u, the top plate is given by

G(u) := {(x, z) ∈ D × R : z = u(x)} .
Recalling that the top plate cannot deformed beyond the surface Σ := D × {-H} of the insulated layer Ω 1 , the deformation u is bounded from below by -H. The region Ω 2 (u) between the top plate and the surface of the insulating layer is defined by

Ω 2 (u) := {(x, z) ∈ D × R : -H < z < u(x)} .
The dielectric permittivity is assumed to be a positive constant σ 2 > 0 in Ω 2 (u) and differs in general from σ 1 (•, -H), so that there is a jump discontinuity of the permittivity across the interface

Σ(u) := {(x, -H) : x ∈ D and u(x) > -H} , separating Ω 1 and Ω 2 (u). Observe that Ω 2 (u) is an open subset of D × (-H, ∞), which is connected when the coincidence set C(u) := {x ∈ D : u(x) = -H} . (1.1)
is empty, while it is disconnected and non-Lipschitz (if u is smooth enough) otherwise. The two situations are depicted in Figure 1 and Figure 2, respectively. Independent of whether or not C(u) is empty, the domain

Ω(u) := Ω 1 ∪ Σ(u) ∪ Ω 2 (u) , is Lipschitz (again if u is smooth enough).
w

Ω 1 Ω 1 Ω 2 (w) Ω 2 (w) D Σ(w) Σ(w) z -H -d -H 0 2L C(w) G(w) FIGURE 2.
Geometry of Ω(u) for a state u = w with non-empty coincidence set (blue).

The total energy E(u) of the MEMS device is given by

E(u) = E m (u) + E c (u) + E e (u) , (1.2)
where -the mechanical energy is

E m (u) := β 2 ∂ 2 x u 2 2 + τ 2 + a 4 ∂ x u 2 2 ∂ x u 2 2 ,
including contributions from bending (β > 0), stretching due to axial tension (τ > 0), and self-stretching due to elongation (a > 0). Here, • 2 denotes the norm in L 2 (D); -the contact energy is

E c (u) := D I [-H,∞) (u) dx , where I [-H,∞) denotes the indicator function of the interval [-H, ∞); -the electrostatic energy is E e (u) := - 1 2 Ω(u) σ|∇ψ u | 2 d(x, z)
with ψ u ∈ H 1 (Ω(u)) denoting the electrostatic potential given as the variational solution to the transmission problem

div(σ∇ψ u ) = 0 in Ω(u) , (1.3a) ψ u = σ∂ z ψ u = 0 on Σ(u) , (1.3b) ψ u = h u on ∂Ω(u) , (1.3c
)

with σ = σ 1 in Ω 1 and σ = σ 2 in Ω 2 (u). In (1.3b),
• denotes the jump across the interface Σ(u), while (1.3c) indicates that ψ u satisfies non-homogeneous Dirichlet boundary conditions prescribed by a given function h u with an explicit dependence upon the deformation u, see (2.3) below. The latter is such that h u ≡ 0 along the bottom plate D × {-H -d} and h u ≡ V along the elastic top plate G(u) with positive potential value V > 0 (see assumption (2.3e) below). The modeling assumption is then that the evolution of u is governed (at least formally) by the gradient flow associated with E, which reads

∂ t u + ∂ u E(u) = 0 , t > 0 , u(0) = u 0 , (1.4)
supplemented with clamped boundary conditions; that is,

u(t) ∈ H 2 D (D) := v ∈ H 2 (D) : v(±L) = ∂ x v(±L) = 0 , t ≥ 0 .
As already noted in [START_REF]Shape derivative of the Dirichlet energy for a transmission problem[END_REF]Section 5], where we studied the existence of minimizers of E, the interpretation of equation (1.4) needs some care for several reasons: First, the contact energy involves a non-smooth convex function and it is rather the notion of subdifferential which is appropriate and requires a suitable functional setting. Specifically, we define

S 0 := v ∈ H 2 D (D) : v > -H in D , S0 := v ∈ H 2 D (D) : v ≥ -H in D , and recall that S0 is a closed convex set in H 2 D (D).
The "derivative" of E c with respect to u is then given by the subdifferential ∂I S0 (u) of the indicator function I S0 of the set S0 . It is a subset of the dual space

H -2 (D) := H 2 D (D) ′ of H 2 D (D) and, for v ∈ S0 , it is given by: ξ ∈ ∂I S0 (v) ⇐⇒ ξ, w -v H 2 D ≤ 0 , w ∈ S0 , where •, • H 2 D
denotes the duality pairing between H -2 (D) and H 2 D (D). Second, the electrostatic energy E e (u) depends on u not only through the integral over Ω(u) but also through the solution ψ u to the transmission problem (1.3). Its differentiability is then a tricky and by no means obvious issue but can be handled with the help of shape optimization tools. It follows from the analysis performed in [START_REF]Shape derivative of the Dirichlet energy for a transmission problem[END_REF] that the functional E e at u ∈ S0 has a directional derivative g(u) given by

g(u)(x) :=      σ 2 2 1 + (∂ x u(x)) 2 ∂ z ψ u,2 (x, u(x)) 2 , x ∈ D \ C(u) , σ 1 (x, -H) 2 2σ 2 ∂ z ψ u,1 (x, -H) 2 , x ∈ C(u) , (1.5)
where ψ u,1 := ψ u | Ω 1 and ψ u,1 := ψ u | Ω 2 (u) . In fact, g(u) corresponds to the electrostatic force acting on the elastic plate. It is worth emphasizing here that the derivation of this result owes much to the book by Henrot & Pierre [START_REF]Shape variation and optimization[END_REF] (and its french version [START_REF] Henrot | Variation et optimisation de formes[END_REF]), which has been a constant source of inspiration in our studies of differentiability properties of the electrostatic energy involved in the modeling of MEMS. In fact, for u ∈ S 0 , the coincidence set C(u) defined in (1.1) is empty and the functional E e is actually Fréchet differentiable at u, a feature which is proved along the lines of [10, Sections 5. out that the functional g(u) involves the traces of ∂ z ψ u on G(u) and C(u) × {-H}, which are well-defined only if ψ u is sufficiently regular. However, since Ω(u) is only a Lipschitz domain while Ω 2 (u) might be even non-Lipschitz when C(u) = ∅, the H 2 (Ω 1 )-regularity of ψ u,1 and the H 2 (Ω 2 (u))-regularity of ψ u,2 are not straightforward and a large part of the analysis performed in [START_REF]Shape derivative of the Dirichlet energy for a transmission problem[END_REF] is devoted to this regularity issue.

Since the computation of the derivative of the mechanical energy with respect to u is classical, collecting the outcome of the above discussion yields the following parabolic variational inequality for u:

∂ t u + β∂ 4 x u -(τ + a ∂ x u 2 2 )∂ 2 x u + ∂I S0 (u) ∋ -g(u) in (0, ∞) × D , (1.6a)
supplemented with the constraint

u(t) ∈ S0 , t ≥ 0 , (1.6b)
and the initial condition

u(0) = u 0 , x ∈ D . (1.6c)
Assuming β > 0, we note that (1.6a) is a fourth-order parabolic variational inequality and the main purpose of this paper is to investigate the existence of weak solutions to (1.6) for a suitable class of boundary data h u occurring in (1.3c), see (2.3) below. In the following, we interpret ∂ 4

x v for v ∈ S0 as an element of H -2 (D) by virtue of

∂ 4 x v, φ H 2 D := D ∂ 2 x v∂ 2 x φ dx , φ ∈ H 2 D (D) .
A definition of a weak solution to (1.6) is then as follows.

Definition 1.1. Let β > 0 and u 0 ∈ S0 . A weak solution to (1.6) is a function u ∈ C([0, ∞), H 1 (D)) ∩ L ∞,loc ([0, ∞), H 2 D (D)) ∩ H 1 loc ([0, ∞), L 2 (D)) satisfying (1.6b), (1.6c), the weak formulation D u(t) -u 0 v dx = - t 0 D g u(s) -(τ + a ∂ x u(s) 2 2 )∂ 2 x u(s) v dxds - t 0 D β∂ 2 x u(s)∂ 2 x v dxds - t 0 ζ(s), v H 2 D (D) ds for any t > 0 and v ∈ H 2 D (D), where ζ := -∂ t u -β∂ 4 x u + τ + a ∂ x u 2 2 ∂ 2 x u -g(u) ∈ L 2,loc ([0, ∞), H -2 (D)) (1.7) with ζ(t) ∈ ∂I S0 u(t) for a.a. t ≥ 0 , (1.8)
and the energy inequality

1 2 t 0 ∂ t u(s) 2 2 ds + E u(t) ≤ E u 0 , t > 0 . (1.9)
The main result of this paper is then the following existence result.

Theorem 1.2. Let β > 0 and τ, a ≥ 0. Suppose that the functions σ and h satisfy, respectively, (2.1) and (2.3) below. Then, given u 0 ∈ S0 , there is a weak solution to (1.6) in the sense of Definition 1.1.

Owing to the variational structure (1.4) of (1.6), the proof of Theorem 1.2 is performed with the help of a time implicit Euler scheme. Given a time step δ > 0, we construct by induction a sequence (u δ n ) n≥0 such that u δ 0 := u 0 and, for n ≥ 0, u δ n+1 is a minimizer of the auxiliary functional

F δ n (v) := 1 2δ v -u δ n 2 + E(v) , v ∈ S0 .
Since F δ n includes a negative contribution from the electrostatic energy E e , we begin the proof by showing that F δ n is bounded below, provided δ ∈ (0, δ 0 ) is sufficiently small, the smallness condition depending only on D, the parameters H, d,β, τ , a, the permittivity σ, and the function h defining the boundary data in (1.3c). The existence of a minimizer u δ n+1 of F δ n on S0 then relies on the lower semicontinuity of the convex part E m + E c of the energy and the properties of the electrostatic energy E e established in [START_REF]Shape derivative of the Dirichlet energy for a transmission problem[END_REF]. As a consequence of u δ n+1 being a minimizer of F δ n on S0 , we further derive a handful of estimates on (u δ n ) n≥1 , which allows us to show that the family (u δ ) δ∈(0,δ 0 ) of piecewise constant functions in time defined by

u δ (t) := ∞ n=0 u δ n 1 [nδ,(n+1)δ) (t) , t ≥ 0 ,
has the right compactness properties, so that its cluster points as δ → 0 are weak solutions to (1.6) in the sense of Definition 1.1. Here again, a key ingredient in the proof is the continuity of the functional g defined in (1.5), which we established in [START_REF]Shape derivative of the Dirichlet energy for a transmission problem[END_REF], see Lemma 3.2 below. Finally, we address the regularity of the distribution ζ ∈ ∂I S0 (u) associated with a weak solution u to (1.6) and given by (1.7)-(1.8). As already mentioned, since β > 0, equation (1.6) is a fourth-order parabolic variational inequality and, as such, the regularity of ζ(t) stemming from (1.8) is that it is a distribution in H -2 (D) for a.e. t > 0. In fact, since the seminal work [START_REF] Ézis | Problèmes unilatéraux[END_REF], regularity for the obstacle problem for the biharmonic parabolic equation has received less attention than the same issue for the obstacle problem for second-order parabolic equations. The only regularity result regarding the obstacle problem for the biharmonic parabolic equation we are aware of is [START_REF] Novaga | Regularity of the obstacle problem for the parabolic biharmonic equation[END_REF], whereas [START_REF] Caffarelli | The obstacle problem for the biharmonic operator[END_REF][START_REF] Frehse | On the regularity of the solution of the biharmonic variational inequality[END_REF][START_REF] Pozzolini | A stability result concerning the obstacle problem for a plate[END_REF][START_REF] Schild | On the coincidence set in biharmonic variational inequalities with thin obstacles[END_REF] are devoted to the elliptic analogue. As in [START_REF] Novaga | Regularity of the obstacle problem for the parabolic biharmonic equation[END_REF], we can prove that -ζ(t) is actually a non-negative bounded Radon measure on D for a.e. t > 0. Let us finally describe the contents of this paper: in the next section, we state the assumptions on the permittivity σ and the boundary data h u in (1.3c). In Section 3, we recall the well-posedness of the transmission problem (1.3) and the regularity of its solution established in [START_REF]Shape derivative of the Dirichlet energy for a transmission problem[END_REF], along with the properties of E e and g from [START_REF]Shape derivative of the Dirichlet energy for a transmission problem[END_REF] which are needed for the analysis performed below. We also show in this section the existence of a minimizer u δ n+1 of the functional F δ n on S0 for sufficiently small values of the time step δ > 0. After this preparation, we are in a position to prove Theorem 1.2 in Section 4 and Corollary 1.3 in Section 5.

ASSUMPTIONS

We provide now the detailed assumptions we put on the permittivity σ and the boundary data h u occurring in the transmission problem (1.3). As already mentioned, the dielectric properties of the device are accounted for by the permittivity σ, which is defined by

σ(x, z) := σ 1 (x, z) for (x, z) ∈ Ω 1 , σ 2 for (x, z) ∈ D × (-H, ∞) ,
where

σ 1 ∈ C 2 Ω 1 with σ 1 > 0 in Ω 1 , σ 2 ∈ (0, ∞) . (2.1)
In particular, there are 0 < σ min < σ max such that

σ min ≤ σ(x, z) ≤ σ max , (x, z) ∈ D × [-H, ∞) . (2.2)
We fix C 2 -functions

h 1 : D × [-H -d, -H] × [-H, ∞) → [0, ∞) (2.3a) and h 2 : D × [-H, ∞) × [-H, ∞) → [0, ∞) (2.3b) satisfying h 1 (x, -H, w) = h 2 (x, -H, w) , (x, w) ∈ D × [-H, ∞) , (2.3c) σ 1 (x, -H)∂ z h 1 (x, -H, w) = σ 2 ∂ z h 2 (x, -H, w) , (x, w) ∈ D × [-H, ∞) .
(2.3d)

Moreover, we assume that

h 1 (x, -H -d, w) = h 2 (x, w, w) -V = 0 , (x, w) ∈ D × [-H, ∞) , (2.3e)
where V > 0, and that there are constants m i > 0, i = 1, 2, 3, such that

|∂ x h 1 (x, z, w)| + |∂ z h 1 (x, z, w)| ≤ m 1 + m 2 w 2 , |∂ w h 1 (x, z, w)| ≤ √ m 3 , (2.3f)
for (x, z, w)

∈ D × [-H -d, -H] × [-H, ∞) and |∂ x h 2 (x, z, w)| + |∂ z h 2 (x, z, w)| ≤ m 1 + m 2 w 2 H + w , |∂ w h 2 (x, z, w)| ≤ m 3 H + w , (2.3g) for (x, z, w) ∈ D × [-H, ∞) × [-H, ∞).
A typical example for a function h satisfying the assumptions (2.3) above was given in [12, Example 5.5] which we recall now.

Example 2.1. Let us consider the situation where σ 1 does not depend on the vertical variable z; that is, σ 1 = σ 1 (x). In that case, we set

h 1 (x, z, w) := V σ 2 (H + z + d) σ 2 d + σ 1 (x)(H + w) , (x, z, w) ∈ D × [-H -d, -H] × [-H, ∞) , and
h 2 (x, z, w) := V σ 2 d + σ 1 (x)(H + z) σ 2 d + σ 1 (x)(H + w) , (x, z, w) ∈ D × [-H, ∞) × [-H, ∞) .
Then assumptions (2.3) are easily checked.

For a given function v ∈ S0 we then define

h v (x, z) := h v,1 (x, z) := h 1 (x, z, v(x)) , (x, z) ∈ Ω 1 , h v,2 (x, z) := h 2 (x, z, v(x)) , (x, z) ∈ D × [-H, ∞) . (2.4)
Let us point out that assumption (2.3c)-(2.3d) guarantee that h v defined in (2.4) satisfies the transmission conditions (1.3b), that is,

h v = σ∂ z h v = 0 on Σ(v) ,
while assumption (2.3e) along with (1.3c) entails that the electrostatic potential ψ v equals zero on the bottom plate D × {-H -d} and equals V along the elastic plate G(u). Assumptions (2.3f)-(2.3g) are required to guarantee the coercivity of the total energy E(v).

Throughout the paper, c and (c i ) i≥1 denote positive constants depending only on D, H, d, β, τ , a, σ, (m i ) 1≤i≤3 , and u 0 . The dependence upon additional parameters will be indicated explicitly.

AUXILIARY RESULTS

We first recall some results that were derived in [START_REF]Shape derivative of the Dirichlet energy for a transmission problem[END_REF] and begin with the well-posedness of (1.3).

Lemma 3.1. [12, Theorem 1.1] Suppose (2.3). For each v ∈ S0 , there is a unique variational solution ψ v ∈ h v + H 1 0 (Ω(v)) to (1.3). Moreover, ψ v,1 := ψ v | Ω 1 ∈ H 2 (Ω 1 ), ψ v,2 := ψ v | Ω 2 (v) ∈ H 2 (Ω 2 (v))
, and ψ v is a strong solution to the transmission problem (1.3)

satisfying σ∂ z ψ v ∈ H 1 (Ω(v)).
The regularity of ψ v stated in Lemma 3.1 guarantees that g(v) defined in (1.5) is meaningful for v ∈ S0 . We collect in the next result some properties of g established in [START_REF]Shape derivative of the Dirichlet energy for a transmission problem[END_REF].

Lemma 3.2. Suppose (2.3). (a) If v ∈ S0 and (v j ) j≥1 ⊂ S0 are such that v j ⇀ v in H 2 (D), then lim j→∞ g(v j ) -g(v) 2 = 0 and lim j→∞ E e (v j ) = E e (v)
. We next turn to differentiability properties of E e . As observed in [START_REF]Shape derivative of the Dirichlet energy for a transmission problem[END_REF], E e need not be Fréchet differentiable for all u ∈ S0 but it has always directional derivatives.

lim s→0 + 1 s E e (v + s(w -v)) -E e (v) = D g(v)(w -v) dx .
We also derive a lower bound on E.

Lemma 3.4. Suppose (2.3).

There is a constant c 1 > 0 such that

E(v) ≥ β 4 ∂ 2 x v 2 2 -c 1 1 + v 2 2 , v ∈ S0 .
Proof. For the sake of completeness, we recall the proof performed in [START_REF]Shape derivative of the Dirichlet energy for a transmission problem[END_REF]Theorem 5.1]. Since ψ v is a variational solution to (1.3), it follows from (2.2), (2.3), and Young's inequality that

-E e (v) = 1 2 Ω(v) σ|∇ψ v | 2 d(x, z) ≤ 1 2 Ω(v) σ|∇h v | 2 d(x, z) ≤ Ω(v) σ (∂ x h(x, z, v(x))) 2 + (∂ w h(x, z, v(x))) 2 (∂ x v(x)) 2 d(x, z) + 1 2 Ω(v) σ (∂ z h(x, z, v(x))) 2 d(x, z) ≤ (d + 1)σ max D 3 2 (m 1 + m 2 v(x) 2 ) + m 3 (∂ x v(x)) 2 dx .
We next use Poincaré's inequality

w 2 ≤ |D| ∂ x w 2 , w ∈ H 1 0 (D) , (3.1)
and the interpolation inequality

∂ x w 2 2 ≤ w 2 ∂ 2 x w 2 , w ∈ H 2 (D) ∩ H 1 0 (D) , (3.2)
to obtain

-E e (v) ≤ c 1 + ∂ x v 2 2 ≤ c 1 + v 2 ∂ 2 x v 2 .
Consequently, Young's inequality and the above upper bound yield

E(v) ≥ β 2 ∂ 2 x v 2 2 -c -c v 2 ∂ 2 x v 2 ≥ β 4 ∂ 2 x v 2 2 -c -c v 2 2 ,
and the proof is complete.

We next provide the basis for the time implicit scheme used later in order to construct a solution to (1.6). Lemma 3.5. Set δ 0 := min{1, (16c 1 ) -1 } > 0. Then, for any δ ∈ (0, δ 0 ) and f ∈ S0 , there is v ∈ S0 such that

- 1 δ (v -f ) -β∂ 4 x v + (τ + a ∂ x v 2 2 )∂ 2 x v -g(v) ∈ ∂I S0 (v) .
Moreover,

1 2δ v -f 2 2 + E(v) ≤ E(f ) .
Proof. The proof relies on the direct method of calculus of variations. Consider δ ∈ (0, δ 0 ) and f ∈ S0 and define

F(v) := 1 2δ v -f 2 2 + E(v) , v ∈ S0 .
Then, by Lemma 3.4 and Young's inequality,

F(v) ≥ 1 2δ v 2 2 2 -f 2 2 + β 4 ∂ 2 x v 2 2 -c 1 1 + v 2 2 ≥ β 4 ∂ 2 x v 2 2 + 1 4δ -c 1 v 2 2 -c 1 - f 2 2 2δ ≥ β 4 ∂ 2 x v 2 2 + 3c 1 v 2 2 -c 1 - f 2 2 2δ
for all v ∈ S0 . Thus, F is bounded from below on S0 and there is a minimizing sequence (v j ) j≥1 in S0 satisfying

µ := inf v∈ S0 {F(v)} ≤ F(v j ) ≤ µ + 1 j , j ≥ 1 . (3.3)
Moreover, the previous lower bound on F guarantees that (v j ) j≥1 is bounded in H 2 D (D). Therefore, there is v ∈ H 2 D (D) such that (up to a subsequence)

v j ⇀ v in H 2 (D) , (3.4a) v j -→ v in H 1 0 (D) . ( 3
E m (v) ≤ lim inf j→∞ E m (v j ) .
Consequently,

F(v) ≤ lim inf j→∞ F(v j ) = µ ,
and we conclude that v ∈ S0 is a minimizer of F on S0 . This property, in turn, guarantees that, for w ∈ S 0 ,

0 ≤ lim inf s→0 + 1 s F(v + s(w -v)) -F(v) .
It then follows from Lemma 3.3 that

0 ≤ D v -f δ (w -v) + β∂ 2 x v ∂ 2 x (w -v) -τ + a ∂ x v 2 2 ∂ x v∂ x (w -v) dx + D g(v)(w -v) dx
for all w ∈ S 0 . Since S 0 is dense in S0 , this inequality also holds for any w ∈ S0 . Therefore,

- 1 δ (v -f ) -β∂ 4 x v + (τ + a ∂ x v 2 2 )∂ 2 x v -g(v) ∈ ∂I S0 (v) .
Finally, since f ∈ S0 , we have F(v) ≤ F(f ), which completes the proof.

PROOF OF THEOREM 1.2

Fix u 0 ∈ S0 . For δ ∈ (0, δ 0 ), we set u δ 0 := u 0 and, using Lemma 3.5, we construct by induction a sequence

(u δ n , ζ δ n ) n≥1 in S0 × H -2 (D) such that ζ δ n+1 := -A δ n+1 -β∂ 4 x u δ n+1 + (τ + a ∂ x u δ n+1 2 2 )∂ 2 x u δ n+1 -g(u δ n+1 ) ∈ ∂I S0 (u δ n+1 ) , (4.1) where A δ n+1 := 1 δ (u δ n+1 -u δ n ) ∈ H 2 D (D) , and
1 2δ u δ n+1 -u δ n 2 2 + E(u δ n+1 ) ≤ E(u δ n ) (4.2)
for n ≥ 0. Let us first note that (4.2) implies

1 2δ n j=0 u δ j+1 -u δ j 2 2 + E u δ n+1 ≤ E(u 0 ) , n ≥ 0 . (4.3)
A first consequence of (4.3) is an L 2 -estimate on (u δ n ) n≥1 , which is adapted from [3, Lemma 3.2.2]. More precisely, it follows from Hölder's and Young's inequalities that, for n ≥ 0,

u δ n+1 2 2 -u 0 2 2 = n j=0 u δ j+1 2 2 -u δ j 2 2 = n j=0 D u δ j+1 -u δ j u δ j+1 + u δ j dx ≤ n j=0 u δ j+1 -u δ j 2 u δ j+1 + u δ j 2 ≤ 1 4c 1 δ n j=0 u δ j+1 -u δ j 2 2 + c 1 δ n j=0 u δ j+1 + u δ j 2 2 .
We then infer from Lemma 3.4 and (4.3) that

u δ n+1 2 2 -u 0 2 2 ≤ E(u 0 ) -E(u δ n+1 ) 2c 1 + 2c 1 δ n j=0 u δ j+1 2 2 + u δ j 2 2 ≤ E(u 0 ) + c 1 + c 1 u δ n+1 2 2 2c 1 + 2c 1 δ u 0 2 2 + 4c 1 δ n+1 j=1 u δ j 2 2 .
Hence,

u δ n+1 2 2 ≤ (2 + 4c 1 δ) u 0 2 2 + 1 + E(u 0 ) c 1 + 8c 1 δ n+1 j=1 u δ j 2 2 , n ≥ 0 .
Since δ ∈ (0, δ 0 ), we have 8c 1 δ < 1/2 < 1 and we are thus in a position to apply a discrete version of Gronwall's lemma, see [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF]Lemma 3.2.4], to conclude that

u δ n+1 2 2 ≤ 6 u 0 2 2 + 2 + 2E(u 0 ) c 1 e 16c 1 (n+1)δ , n ≥ 0 . (4.4)
We next use again Lemma 3.4 and (4.3), along with (4.4), to obtain that, for n ≥ 0,

1 2δ n j=0 u δ j+1 -u δ j 2 2 + β 4 ∂ 2 x u δ n+1 2 2 ≤ 1 2δ n j=0 u δ j+1 -u δ j 2 2 + E u δ n+1 + c 1 1 + u δ n+1 2 2 ≤ E(u 0 ) + c 1 1 + u δ n+1 2 2
≤ c 2 e 16c 1 nδ .

Owing to the functional inequalities (3.1) and (3.2), we end up with

1 2δ n j=0 u δ j+1 -u δ j 2 2 + u δ n+1 2 H 2 ≤ c 3 e 16c 1 nδ , n ≥ 0 . (4.5)
We next define the functions

u δ , A δ : [0, ∞) → H 2 D (D) and ζ δ : [0, ∞) → H -2 (D) by u δ (t) := n≥0 u δ n 1 [nδ,(n+1)δ) (t) , t ≥ 0 , A δ (t) := n≥1 A δ n 1 [nδ,(n+1)δ) (t) , t ≥ 0 , and ζ δ (t) := n≥1 ζ δ n 1 [nδ,(n+1)δ) (t) , t ≥ 0 ,
respectively.

Lemma 4.1. There are a sequence δ ℓ → 0 and

u ∈ C([0, ∞), H 1 (D)) ∩ L ∞ ((0, ∞), H 2 D (D)) ∩ H 1 loc ([0, ∞), L 2 (D)) with u(0) = u 0 such that u(t) ∈ S0 , t ≥ 0 , (4.6
) and, for all t > 0,

u δ ℓ (t) → u(t) in H 1 (D) , (4.7a)
u δ ℓ ⇀ u in L 2 ((0, t), H 2 D (D)) , (4.7b)
g u δ ℓ → g(u) in L 2 ((0, t), L 2 (D)) , (4.7c)
A δ ℓ ⇀ ∂ t u in L 2 ((0, t), L 2 (D)) . (4.7d)
In particular,

ζ δ ℓ ⇀ ζ in L 2 ((0, t), H -2 (D))
for any t > 0 and

ζ := -∂ t u -β∂ 4 x u + τ + a ∂ x u 2 2 ∂ 2 x u -g(u) ∈ L 2,loc ([0, ∞), H -2 (D)) . (4.8)
Proof. Given 0 ≤ t 1 < t 2 there are integers

n 1 ≤ n 2 such that t i ∈ [n i δ, (n i + 1)δ), i = 1, 2.
Either n 1 = n 2 and u δ (t 2 ) = u δ (t 1 ). Or n 1 < n 2 and we infer from (4.5) and Hölder's inequality that

u δ (t 2 ) -u δ (t 1 ) 2 = u δ n 2 -u δ n 1 2 ≤ n 2 -1 j=n 1 u δ j+1 -u δ j 2 ≤ √ n 2 -n 1   n 2 -1 j=n 1 u δ j+1 -u δ j 2 2   1/2 ≤ √ 2c 3 n 2 δ -n 1 δ e 8c 1 n 2 δ . Thus, u δ (t 2 ) -u δ (t 1 ) 2 ≤ √ 2c 3 t 2 -t 1 + δ e 8c 1 t 2 , 0 ≤ t 1 < t 2 . (4.9)
Moreover, for t > 0, there is an integer n ≥ 0 such that t ∈ [nδ, (n + 1)δ) and, either n = 0 and u δ (t) = u 0 , or, again by (4.5), u δ (t) Since H 2 (D) embeds compactly in L 2 (D), we infer from (4.9) and (4.10) that we may apply a variant of the Arzelà-Ascoli theorem, see [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF]Proposition 3.3.1], and a diagonal argument to obtain the existence of a sequence (δ ℓ ) ℓ≥1 , δ ℓ → 0, and

u ∈ C([0, ∞), L 2 (D)) ∩ L ∞,loc ([0, ∞), H 2 D (D)) such that u δ ℓ (t) → u(t) in L 2 (D) , t ≥ 0 , (4.11)
and

u δ ℓ ⇀ u in L 2 ((0, T ), H 2 D (D)) , T > 0 .
We have thus proved (4.7b). Next, an interpolation argument, together with (4.10) and (4.11), yields (4.7a) and the stated time continuity of u in H 1 (D). Furthermore, combining (4.7a), (4.10), and Lemma 3.2 allows one to apply Lebesgue's theorem to deduce (4.7c). Also, since u δ (0) = u 0 and u δ (t) ∈ S0 for t ≥ 0 and δ ∈ (0, δ 0 ), we readily deduce from (4.11) that u(0) = u 0 and u(t) ∈ S0 for t ≥ 0.

Next, for t > 0, there is n ℓ ≥ 0 such that t ∈ [n ℓ δ ℓ , (n ℓ + 1)δ ℓ ) and it follows from (4.5) that

t 0 A δ ℓ (s) 2 2 ds ≤ (n ℓ +1)δ ℓ 0 A δ ℓ (s) 2 2 ds = 1 δ ℓ n ℓ j=0 u δ j+1 -u δ j 2 2 ≤ 2c 3 e 16c 1 t . (4.12)
Since

A δ (t, x) = u δ (t, x) -u δ (t -δ, x) δ , (t, x) ∈ (δ, ∞) × D ,
and A δ (t, x) = 0 for (t, x) ∈ (0, δ)×D, the sequence (A δ ℓ ) ℓ≥1 converges to ∂ t u in D ′ ((0, ∞)×D) as ℓ → ∞, so that the just established boundedness of (A δ ℓ ) ℓ≥1 in L 2,loc ([0, ∞), L 2 (D)) implies that ∂ t u ∈ L 2,loc ([0, ∞), L 2 (D)) and the convergence (4.7d) (up to a subsequence). The stated convergence of ζ δ ℓ ℓ≥1 and the regularity of ζ are straightforward consequences of the regularity of u and (4.7).

We next prove the energy inequality (1.9). Lemma 4.2. For t > 0,

1 2 t 0 ∂ t u(s) 2 2 ds + E(u(t)) ≤ E(u 0 ) .
Proof. Given t > 0 and ℓ ≥ 0, we pick again the integer n ℓ such that t ∈ [n ℓ δ ℓ , (n ℓ + 1)δ ℓ ). Then, by (4.3),

1 2δ ℓ n ℓ j=0 u δ j+1 -u δ j 2 2 + E(u δ ℓ (t)) ≤ E(u 0 ) . ( 4
-A δ n + (τ + a ∂ x u δ n 2 2 )∂ 2 x u δ n -g(u δ n ) u δ n -v dx -β D ∂ 2 x u δ n ∂ 2 x u δ n -v) dx (4.17) for n ≥ 1, δ ∈ (0, δ 0 ), and v ∈ S0 . Now, consider a non-negative function φ ∈ C c ([0, ∞)) and w ∈ L 2,loc ([0, ∞), H 2 (D)
) such that w(t) ∈ S0 for a.e. t ∈ (0, ∞). Then, for δ small enough, supp φ ⊂ (δ, ∞) and, by (4.17),

∞ 0 φ(s) D -A δ (s) + τ + a ∂ x u δ (s) 2 2 ∂ 2 x u δ (s) -g(u δ (s)) u δ (s) -w(s) dxds - ∞ 0 φ(s) D β∂ 2 x u δ (s)∂ 2 x u δ (s) -w(s)) dxds = ∞ n=1 (n+1)δ nδ φ(s) D -A δ n + τ + a ∂ x u δ n 2 2 ∂ 2 x u δ n -g(u δ n ) u δ n -w(s) dxds - ∞ n=1 (n+1)δ nδ φ(s) D β∂ 2 x u δ n ∂ 2 x u δ n -w(s)) dxds ≥ 0 . (4.18)
On the one hand, it follows from (4.7) that

lim ℓ→∞ ∞ 0 φ(s) D -A δ ℓ (s) + (τ + a ∂ x u δ ℓ (s) 2 2 )∂ 2 x u δ ℓ (s) -g(u δ ℓ (s)) u δ ℓ (s) -w(s) dxds = ∞ 0 φ(s) D -∂ t u(s) + (τ + a ∂ x u(s) 2 2 )∂ 2 x u(s) -g(u(s)) u(s) -w(s) dxds .
On the other hand, we infer from (4.7b) and the non-negativity of φ that

lim ℓ→∞ ∞ 0 φ(s) D β∂ 2 x u δ ℓ (s)∂ 2 x w(s) dxds = ∞ 0 φ(s) D β∂ 2 x u(s)∂ 2 x w(s) dxds and lim inf ℓ→∞ ∞ 0 φ(s) D β|∂ 2 x u δ ℓ (s)| 2 dxds ≥ ∞ 0 φ(s) D β|∂ 2 x u(s)| 2 dxds .
Collecting the above identities, taking δ = δ ℓ in (4.18), and letting ℓ → ∞, we conclude that

∞ 0 φ(s) D -∂ t u(s) + (τ + a ∂ x u(s) 2 2 )∂ 2 x u(s) -g(u(s)) u(s) -w(s) dxds - ∞ 0 φ(s) D β∂ 2 x u(s)∂ 2 x u(s) -w(s)) dxds ≥ 0 . That is, recalling the definition (4.8) of ζ, ∞ 0 φ(s) ζ(s), u(s) -w(s) H 2 D ds ≥ 0 for any w ∈ L 2,loc ([0, ∞), H 2 D (D)
) satisfying w(t) ∈ S0 for a.e. t ∈ (0, ∞) and any nonnegative φ ∈ C c ([0, ∞)). In particular, for all v ∈ S0 and non-negative φ ∈ C c ([0, ∞)), the choice w(t) ≡ v, t > 0, in the above inequality gives

∞ 0 φ(s) ζ(s), u(s) -v H 2 D ds ≥ 0 , which implies, since S0 is separable, that ζ(t) ∈ ∂I S0 u(t) for a.a. t ≥ 0 .

Corollary 1 . 3 .

 13 Let the assumptions of Theorem 1.2 be satisfied. If u is a weak solution to (1.6) in the sense of Definition 1.1, then u ∈ L 2,loc ([0, ∞), H s (D)) for s ∈ (2, 7/2) and -ζ ∈ L 2,loc ([0, ∞), M + (D)), where M + (D) is the positive cone of the space M(D) = C 0 (D) ′ of bounded Radon measures on D.

( b )

 b The mapping g : S0 → L 2 (D) is continuous and bounded on bounded sets, the set S0 being endowed with the topology of H 2 (D). Proof. Since weak convergence in H 2 (D) implies boundedness in H 2 (D) and strong convergence in H 1 (D), part (a) is shown in [12, Proposition 3.17 & Corollary 3.12]. As for part (b), the continuity of g : S0 → L 2 (D) follows from [12, Theorem 1.4], while the boundedness of g on bounded sets is a consequence of [12, Corollary 3.14 & Lemma 3.16] and the continuity of the trace from H 1 (Ω 1 ) to L p (D × {-H}) for all p ∈ [1, ∞).

Lemma 3 . 3 .

 33 [START_REF]Shape derivative of the Dirichlet energy for a transmission problem[END_REF] Proposition 5.6] Suppose (2.3). Let v ∈ S0 and w ∈ S 0 . Then

  .4b) Clearly, (3.4a) ensures that v ∈ S0 since S0 is closed and convex in H 2 D (D). It then follows from (3.4a) and Lemma 3.2 (a) that E e (v) = lim j→∞ E e (v j ) , while (3.4) and the weak lower semicontinuity of the L 2 -norm readily imply that

Finally 2 D 2 ). 3 ) 1 . 2 ≥ 2 ∈ 2 D=. 2 D. 5 ) 2 . 2 D 2 D 2 D 2 D

 22312222522222 , since ∂ t u ∈ L 2,loc ([0, ∞), L 2 (D)), it follows from the definition (4.8) of ζ that u solves(1.6) in the sense of Definition 1.1, and the proof of Theorem 1.2 is complete.5. PROOF OF COROLLARY 1.3We finally derive the additional features enjoyed by ζ as stated in Corollary 1.3.Proof of Corollary 1.3. Let u be a weak solution to(1.6) in the sense of Definition 1.1 and define ζ by (1.7). We introduce the setZ := {t ∈ (0, ∞) : ζ(t) ∈ ∂I S0 u(t) } ,(5.1)and observe that |Z| = 0 since ζ satisfies (1.8). Moreover, since u(t) + v belongs to S0 for any non-negative v ∈ C ∞ c (D), it readily follows from (5.1) that -ζ(t), v H That is, for t ∈ Z, -ζ(t) is a non-negative distribution on D and thus a non-negative Radon measure, see, e.g., [7, Proposition 6.6].Next, let T > 0. According to the regularity of u,K T := sup t∈[0,T ] { u(t) H 2 } < ∞ . (5Step Let t ∈ [0, T ] and x ∈ [0, L]. Since u(t) ∈ H 2 D (D), it follows from (5.3) that u(t, x) = u(t, L) + (x -L)∂ x u(t, L) + x)|∂ 2 x u(t, y)| dy ≥ -(Lx) -K T (Lx) 3/2 .Hence, for x ∈ [x T , L] withx T := max L -[L/2, L) ,we deduce that u(t, x) ≥ -H/2. Using the same argument for x ∈ [-L, 0], we end up withu(t, x) ≥ -H 2 , x ∈ [-L, -x T ] ∪ [x T , L] , t ∈ [0, T ] . (5.4) Now, consider t ∈ [0, T ] ∩ Z and v ∈ C ∞ c (D) such that supp v ⊂ [-L, -x T ] ∪ [x T , L].For θ ∈ (0, 1) small enough, we infer from (5.4) that u(t) ± θv belongs to S0 , so that (1.8) entails0 ≤ ζ(t), u(t)u(t) ∓ θv H ∓θ ζ(t), v H 2 D Since θ is positive, we conclude that ζ(t), v H = 0. Consequently, supp ζ(t) ⊂ [-x T , x T ] , t ∈ [0, T ] ∩ Z . (5Step We now fix a non-negative functionχ T ∈ C ∞ c (D) such that χ T ≡ 1 on [-x T , x T ]. Then, for v ∈ C ∞ c (D) and t ∈ [0, T ]∩ Z, the function u(t)+ χ T ( v ∞ ± v)belongs to S0 and it follows from (5.2) and (5.5) that0 ≤ -ζ(t), χ T ( v ∞ ± v) H = v ∞ -ζ(t), χ T H 2 D ± -ζ(t), v H 2 D , the identity ζ(t), χ T v H = ζ(t), v Hbeing guaranteed by(5.5) and the properties of χ T . Thus,-ζ(t), v H ≤ v ∞ -ζ(t), χ T H 2 D ,and the density ofC ∞ c (D) in C 0 (D) and the already established non-negativity of -ζ ensure that -ζ(t) belongs to M + (D) with ζ(t) M ≤ -ζ(t), χ T H 2 D , t ∈ [0, T ] ∩ Z . (5.6) Since ζ ∈ L 2 ((0, T ), H -2 (D)), an immediate consequence of (5.6) is that ζ ∈ L 2 ((0, T ), M(D)).Finally, let s ∈ (2, 7/2). According to [1, Lemma 4.1 (iii)], M(D) is continuously embedded in H s-4 D (D), so that the just established regularity of ζ implies that ζ ∈ L 2 ((0, T ), H s-4 D (D)). Together with the regularity of u, ∂ t u, and g(u), this property and (1.7) ensure that ∂ 4 x u belongs to L 2 ((0, T ), H s-4 D (D)). Consequently, it follows from elliptic regularity theory that u belongs to L 2 ((0, T ), H s D (D)). This completes the proof of Corollary 1.3. Remark 5.1. Since u ∈ C([0, ∞) × D) by Theorem 1.2 and since H 1 (D) embeds continuously in C( D), it easily follows from (1.8) (by an argument similar to that leading to (5.5)) that supp ζ(t) ⊂ C(u(t)) for a.e. t ∈ (0, ∞) , the coincidence set C(u(t)) being defined in (1.1).

  2 H 2 ≤ c 3 e 16c 1 t . Consequently, u δ (t) H 2 ≤ c 4 e 16c 1 t , t ≥ 0 . (4.10)

  Since (u δ ℓ (t)) ℓ≥0 is bounded in H 2 (D) according to (4.10), we may extract a further subsequence (not relabeled), possibly depending on t, such that u δ ℓ (t) ℓ≥0 converges to u(t) weakly in H 2 (D) and strongly in H 1 (D). Hence Proof of Theorem 1.2. To finish off the proof, we are left with showing that u solves the variational inequality (1.8). To this end, we recall from (4.1) that

					E m u(t) ≤ lim inf ℓ→∞	E m u δ ℓ (t) ,
	which gives, together with (4.14),
					E u(t) ≤ lim inf ℓ→∞	E u δ ℓ (t) .	(4.15)
	Moreover, due to (4.7d) and (4.12), we have
		1 2	0	t	∂ t u(s) 2 2 ds ≤ lim inf ℓ→∞	1 2δ ℓ	n ℓ j=0	u δ ℓ j+1 -u δ ℓ j	2 2 .	(4.16)
	Gathering (4.13), (4.15), and (4.16) completes the proof.
	0 ≤	D
										.13)
	Owing to Lemma 3.2 (a), (4.7a), and (4.7b), we have
					lim ℓ→∞	E e u δ ℓ (t) = E e u(t) .	(4.14)


Olga Rodionova
email: rodionova@recherche.enac.fr

Daniel Delahaye
email: delahaye@recherche.enac.fr

Mohamed Sbihi
email: sbihi@recherche.enac.fr

Marcel Mongeau
email: mongeau@recherche.enac.fr

TRAJECTORY PREDICTION IN NORTH ATLANTIC OCEANIC AIRSPACE BY WIND NETWORKING

Aircraft in the North Atlantic oceanic airspace are subjected to strong winds caused by the Jet Stream. The flight progress prediction is based on meteorological data forecast. The roughness in wind forecast often results in large differences between the predicted and the actual cruising times. This increases the uncertainty in conflict detection. All aircraft are able to perform instant meteorological measurements. Implementing new technologies enables exchanging these measured data directly between aircraft. Thus, an aircraft can obtain information about winds on its route from preceding aircraft and adjust its time predictions. This process is called Wind Networking. As the data obtained with wind networking is much more accurate than the initial forecast, the adjusted predictions are much closer to the reality.

In this work, we perform simulations of the wind networking process. We compare the results of applying this approach and the results obtained using the meteorological forecast in terms of cruising time prediction and conflict prediction. We observe a remarkable improvement in the predictions by the wind networking and we quantify the benefits of using such an approach.

Introduction

The North Atlantic oceanic airspace (NAT) accommodates air traffic between two densely populated areas -North America and Europe. Aircraft crossing NAT are subjected to very strong winds caused by the presence of Jet Streams (JSs). JSs are fast narrow predominantly west-east air currents mainly located in the upper troposphere that are caused by a combination of the earth rotation and atmospheric heating (Figure 1). Two types of JSs are distinguished in NAT: subtropical JSs raised by the westerly acceleration of poleward moving air in the upper tropical circulation and eddy-driven JSs provoked by the momentum and heat forcing from transient midlatitude eddies [START_REF] Woollings | Variability of the North Atlantic eddy-driven jet stream[END_REF][START_REF] Gerber | On the zonal structure of the North Atlantic oscillation and annular modes[END_REF]. The speed of the JS is typically 100 kts (nautical miles per hour) but can reach 200 kts.

Figure 1. Jet Stream in NAT 1

Air traffic in NAT mainly contributes to two major flows: the westbound flow departing from Europe in the morning, and the eastbound flow departing from North America in the evening. Because of the time zone differences and passenger demands, these flows are separated in time. In the free flight conditions, the eastbound flights would try to exploit the JS in order to benefit from strong tailwinds, while the westbound flights would avoid the JS in order to minimize headwind [START_REF] Irvine | Characterizing North Atlantic weather patterns for climate-optimal aircraft routing[END_REF]. As a result the two main flows are found to be separated in space as well; thus, they can be treated independently.

There are more than 500 flights crossing the NAT in each direction daily. As their wind-optimal routes are very close to each other, the NAT is highly congested at peak hours. Standard radar-based surveillance is not available for most part of this vast airspace. In order to increase efficiency of the Air Traffic Control (ATC) [START_REF] Wu | Performance metrics for oceanic air traffic management[END_REF] and to perform conflictfree flight progress, the concept of parallel flight tracks known as the Organized Track System (OTS) was introduced in NAT [START_REF]Guidance concerning air navigation in and above the North Atlantic MNPS airspace[END_REF]. Two independent OTS for westbound and eastbound flights are constructed daily based on the minimal-time aircraft routes depending on the position of anticyclones and jet winds (Figure 2). Mainly as a result, the eastbound tracks are located more northerly than the westbound tracks. The eastbound (night-time) OTS is valid from 0100 UTC (Coordinated Universal Time, hours and minutes of the Greenwich mean time) to 0800 UTC, while the hours of validity of the westbound (daytime) OTS are from 1130 UTC to 1900 UTC.

Figure 2. An instance of the OTS 2

The OTS consists of several (from 5 to 8) quasiparallel tracks constructed so that the lateral separation of 60 NM (nautical miles) is automatically maintained for aircraft flying on adjacent tracks (Figure 3). Each track represents a sequence of great circles joining successive significant waypoints (WPs). The WPs include named points and latitude crossings of all oceanic ten-degree meridians at integer degrees of latitude. Most adjacent tracks are separated by 1º of latitude at each WP. In the vertical direction, each track consists of several flight levels (FLs) between FL290 and FL410 inclusive (Figure 3). The distance between adjacent FLs is 1000 feet, in order to ensure the minimal required vertical separation.

Currently, about half of NAT flights use the OTS [START_REF]Guidance concerning air navigation in and above the North Atlantic MNPS airspace[END_REF]. An aircraft intended to use the OTS should enter a predefined track at its entry WP at a predefined FL and follow this track at the same FL until any modification of its planned trajectory is required and authorized. An aircraft can change its speed (Mach number) or/and its FL only at WPs and only after acquiring the corresponding Clearance (permission) from the Oceanic Area Control Center (OAC) [START_REF]Guidance concerning air navigation in and above the North Atlantic MNPS airspace[END_REF]. An aircraft is also authorized to demand the re-routing from its current track to an adjacent one, also only at WPs, but this maneuver is very rarely applied in NAT because of the high traffic density and the large longitudinal separation standards [START_REF] Louyot | ASEP-ITM simulations from traffic data, DSNA[END_REF][START_REF] Rodionova | Optimization of aicraft trajectories in North Atlantic oceanic airspace[END_REF].

Figure 3. NAT traffic model 3

The longitudinal separation (between aircraft flying on the same track) in NAT is assessed in terms of differences between actual and estimated times of arrival at common WPs. It is expressed in clock minutes. This time-based separation is defined so as to guarantee the separation minima whenever en route. According to the current separation standards (CSS) in NAT [START_REF]Guidance concerning air navigation in and above the North Atlantic MNPS airspace[END_REF], two consecutive aircraft following the same track should be separated at least 10 minutes apart (in some cases this norm can be reduced down to 5 minutes [START_REF][END_REF]). In the case of rerouting of a particular aircraft onto an adjacent track, an increased longitudinal separation of 15 minutes is imposed [START_REF] Louyot | ASEP-ITM simulations from traffic data, DSNA[END_REF] because of the difference in expected meteorological conditions between the tracks (Figure 4).

Implementing new technologies [START_REF] Mchale | NextGen air traffic control avionics are moving from concept to cockpit[END_REF] (such as the Automatic Dependent Surveillance-Broadcast (ADS-B) technology [START_REF] Mchale | ADS-B in enhances ATM and pilot situational awareness[END_REF][START_REF] Cho | The improvement of aircraft position information with the unscented Kalman filter[END_REF]) permits to apply the airbornedelegated separation to aircraft [START_REF] Peng | Study on the model for horizontal escape maneuvers in TCAS[END_REF][START_REF] Prandini | Toward air traffic complexity assessment in new generation air traffic management systems[END_REF]. This in turns is expected to lead to a significant decrease of the CSS [START_REF] Massoglia | The use of satellite technology for oceanic air traffic control[END_REF]. Indeed, the reduced separation standards (RSS) should reach 2 minutes for consecutive aircraft on the same track, and 3 minutes for re-routing aircraft [START_REF] Louyot | ASEP-ITM simulations from traffic data, DSNA[END_REF][START_REF] Rodionova | Optimization of aicraft trajectories in North Atlantic oceanic airspace[END_REF]. Using the RSS raises significantly the NAT traffic accommodation capacity and the efficiency of individual flights by allowing the aircraft to follow more optimal routes [START_REF] Williams | Benefits assessment of reduced separations in North Atlantic Organized Track System[END_REF][START_REF] Rodionova | North Atlantic Aircraft Trajectory Optimization[END_REF]. At the same time, it needs precise aircraft position determination and thus, more precise trajectory prediction.

Figure 4. Longitudinal separation standards

Currently, the flight progress prediction relies on the meteorological forecast (MF) data (especially, wind maps). Typically MF data is not sufficiently precise [START_REF] Hollingsworth | The statistical structure of shortrange forecast errors as determined from radiosonde data[END_REF][START_REF] Schwartz | Accuracy of RUC-1 and RUC-2 wind and aircraft trajectory forecasts by comparison with ACARS observations[END_REF][START_REF] Cardinali | Jet stream analysis and forecast errors using GADS aircraft observations in the DAO, ECMWF, and NCEP models[END_REF][START_REF] Mondoloni | Aircraft trajectory prediction errors: Including a summary of error sources and data[END_REF] and yields substantial errors in aircraft trajectory predictions [START_REF]Navigational errors on the North Atlantic, the Global Voice of Pilots[END_REF][START_REF] Lee | Performance of trajectory models with wind uncertainty[END_REF]. This in turn increases the uncertainty in conflict detection and the difficulty of tactical conflict resolution. At the same time, all aircraft are supplied with equipment capable to perform instant meteorological measurements (temperature, pressure, wind speed) during the flight [START_REF] Wood | Aircraft observations of boundary layer structure[END_REF][START_REF] Khelif | Improved wind measurements on research aircraft[END_REF][START_REF] Moninger | Automated meteorological reports from commercial aircraft[END_REF]. Even though such measurements are always carried out with some inaccuracy [START_REF] Benjamin | Accuracy of ACARS wind and temperature observations determined by collocation[END_REF][START_REF] Drue | Aircraft type-specific errors in AMDAR weather reports from commercial aircraft[END_REF][START_REF] Moninger | Evaluation of regional aircraft observations using TAMDAR[END_REF], having the possibility to use them in wind forecast adjustment can significantly improve flight predictions [START_REF] Cole | An assessment of the 60km Rapid Update Cycle (RUC) with near real time aircraft reports[END_REF][START_REF] Cole | Improving RUC-1 wind estimates by incorporating near-realtime aircraft reports[END_REF].

New surveillance and broadcast services (ADS-B, Flight Information Services-Broadcast (FIS-B) [START_REF] Zimmerman | ADS-B 101: what it is and why you should care, Air Facts, the journal for personal air travel -by pilots[END_REF][START_REF] Gilbert | Surveillance and broadcast services -an effective nationwide solution[END_REF]) allow aircraft to exchange the measured meteorological data with the ground ATC and with each other directly en-route (Figure 5). Thus, an aircraft can obtain recent and more precise information about winds on its route from preceding aircraft, and can update its time predictions based on these data. We call this process: Wind Networking (WN).

The current work is devoted to evolve and quantify the benefits of using WN approach in aircraft trajectory prediction and conflict detection. First we describe the modeling of the NAT air traffic situation. Next we discuss more precisely the WN approach. Finally, the simulation results are presented followed by the comparison of the two prediction methods: MF and WN.

NAT air traffic models

This section describes the models we propose for the NAT air traffic simulations: wind model, OTS model and flight model. Further modeling of the real flight progress and its prediction using the MF is introduced. In this work, we focus on simulating the NAT eastbound traffic. Westbound traffic can be treated in the same manner.

Wind model

The quality of aircraft trajectory prediction depends first of all on the quality of wind prediction. To evaluate the prediction accuracy, one needs to compare the predicted wind values with the real ones. The difficulty is that the real wind is actually unknown. The forecast wind can be obtained from meteorological databases. Thus, a large amount of works devoted to modeling wind for aviation focus on modeling the deviations of the predicted wind from the real one [START_REF] Mondoloni | Aircraft trajectory prediction errors: Including a summary of error sources and data[END_REF][START_REF] Glover | A stochastic hybrid model for air traffic control simulation[END_REF][START_REF] Lymperopoulos | Sequential Monte Carlo methods for multi-aircraft trajectory prediction in air traffic management[END_REF]. These deviations are generally given under statistical form [START_REF] Hollingsworth | The statistical structure of shortrange forecast errors as determined from radiosonde data[END_REF][START_REF] Tenenbaum | Determining biases in Hadley circulation reanalyses using independent aircraft observations[END_REF] that requires adjusting several parameters. As to our knowledge there was no precise study similar to [START_REF] Cole | An assessment of the 60km Rapid Update Cycle (RUC) with near real time aircraft reports[END_REF] determining the statistical parameters for NAT, we propose another approach to model both real and forecast winds, and their deviation.

We rely on the GRIB data files created by the meteorological centers (NOAA, Meteoblue, FNMOC), which are available online [START_REF]Weather data visualization[END_REF]. A GRIB file consists of several grids coded in a special format [START_REF] Dey | GRIB (Edition 1) The WMO format for the storage of weather product information and the exchange of weather product messages in gridded binary form as used by NCEP central operations[END_REF]. Each such grid represents a regular grid of geographical points with maximal resolution of 0.5º, where various MF data is defined for each point. We study the area from 35ºN to 70ºN of latitude and from 0ºW to 75ºW of longitude, which covers the NAT. We extract only the u (west-east) and v (north-south) wind components for each point. Each grid is defined for a specific time for which the MF is valid (with 3hour time step) and for a specific altitude. We concentrate on the time period from 0000 UTC to 0900 UTC that includes the eastbound OTS validity period, and on the altitudes 300hPa and 200hPa that bound the published OTS FLs. We decode the GRIB files with the open-source java library JGrib [START_REF]JGRIB[END_REF].

In order to simulate one day of traffic, we rely on two GRIB files downloaded with time difference of several (from 18 to 24) hours. We use the data from the first (the oldest) file as forecast winds. The second (the latest) file is used to model the real winds, as its data is more recent and therefore tends to represent better the reality. Figure 6 displays the wind fields obtained from the two GRIB files available for 10 December 2013 at the altitude 200 hPa. The left column corresponds to the forecast wind; the right column simulates the assumed real wind. Two time moments: 0000 UTC (upper line) and 0600 UTC (lower line) are considered in this example to illustrate the wind evolution. As it can be seen from Figure 6, the forecast and real wind fields differ slightly.

Figure 6. Forecast and real wind fields

GRIB files only provide information about the winds at discrete grid points. We obtain the continuous wind fields by linear interpolation over grids in 4-dimensional space (time, altitude, latitude and longitude). With current grid resolution the 4D linear interpolation coincides satisfactorily with the reality. Remark that most wind models also rely on such an approach.

OTS model

As the OTS is constructed daily based on the wind fields, it is important to use wind forecast and OTS valid corresponding to the same day in the simulations. The current [START_REF]Current NAT tracks[END_REF] and the archived [START_REF]Enroute Flightplanning Resource & Archive[END_REF] OTS can be found online in free access in the form of NAT messages. Such messages contain a number of OTS tracks defined for a particular day, where each track is represented by a sequence of WPs followed by a sequence of FLs available for the current track. Figure 7 displays the eastbound OTS defined on 10 December 2013 projected on the horizontal plane. grid of WPs, where Ny is the number of OTS tracks, Nx is the number of WPs on each track, and Nz is the number of FLs for each track. The tracks are labeled from 1 to Ny starting from the most northern, the WPs on each track from 1 to Nx starting from the most western, and the FLs from 1 to Nzstarting from the lowest. Thus, the 3D position of an aircraft on track j at WP i at FL k is completely specified by the vector  

k j i , ,
. The geographical coordinates of this point are given by the vector

  k j i j i h , ,  , where j i  is the point longitude, j i  is
the point latitude and k h is the point altitude.

Every pair of consecutive WPs on the same track:  

k j i , ,
and  

k j i , , 1 
are joined with straight lines called links. Similarly, every pair of consecutive WPs on adjacent tracks, northern:  

k j i , , and
  k j i , 1 , 1   ; southern:   k j i , ,
and  

k j i , 1 , 1  
), are joint with links. Thus, each WP (except for those on the outer tracks) has three outgoing links. We call the link intersection points as nodes. Figure 8 displays the horizontal section of our OTS model represented by a grid of nodes joined with links. In addition to WPs that are located on tracks and that belong to the OTS grid structure (red nodes in Figure 8) the node set also contains the points of intersection of the links joining adjacent tracks (blue nodes in Figure 8).

Flight model

We considered several flight sets with various numbers of flights in our simulations. For some tests the flight input parameters were set randomly (respecting both the ATC regulations and our model constraints). Other tests were performed using real flight sets extracted from the FPL (flight plan) messages. Such messages are available for several days of August 2006 from the ASSTAR project [START_REF] Louyot | ASEP-ITM simulations from traffic data, DSNA[END_REF] (we did not have the access to more recent data).

Here is how we simulate the flight progress of N aircraft crossing the NAT within the OTS. For each flight f we define the entry ( f in Tr ) and exit ( f out Tr ) tracks, and the track entry time ( f in t ). If the entry and exit tracks are different, we choose the WPs where the aircraft performs re-routings. Thus, the aircraft route between the entry and exit WPs consists of a sequence of links. We limit our model to authorize only re-routings towards the exit tracks in order to avoid north-south zigzagging.

We assume that while moving along a link, an aircraft maintains constant air speed (Mach number) and constant FL. Thus, true air speed ( f i v ) and FL ( f i FL ) are to be set only at WPs ( Nx i ,..., 1 

). Furthermore, we assume that the speed difference between two consecutive WPs ( i and 1  i ) does not exceed 0.02 Mach. Current practice tends to show that the FLs difference is very rarely above 2 FLs [START_REF][END_REF]. In addition to this, aircraft are allowed only to climb (not to descend), in order to satisfy the optimal kerosene consumption flight profile. Figure 9 displays the vertical section of the OTS grid and shows possible altitude profiles for eastbound flights. In our model, the distance between the real horizontal position of the aircraft at the new FL and the previous WP after climbing (see Figure 9) is neglected, as well as the time required to reach the new FL. This instantaneous-climbing hypothesis is reasonable, as the distance between the WPs (10º of longitude, or about 500 km) is much longer than the vertical distance between FLs (1000 feet, or 0.3048 km). The particularity of our OTS model imposes that the flight interactions can only happen on links or at nodes. Lateral and vertical separation being ensured by the OTS construction, a conflict may only occur when a longitudinal separation constraint is violated for a pair of consecutive flights. A conflict can occur on a link if one aircraft is slower than the other following it on the same track. A conflict is detected on a node (shown with a red circle in Figure 10) when two consecutive aircraft pass this node with a time difference violating the longitudinal time separation minimum. The task of ATC is to detect and eliminate such conflicts (see [START_REF] Rodionova | North Atlantic Aircraft Trajectory Optimization[END_REF]). In the current paper we focus only on conflict detection and on the comparison of prediction methods.

Flight progress simulation and prediction

To simulate flight progress within NAT, we need to compute for each aircraft the times of passing the route points. These times are computed based on the aircraft ground speed. While the aircraft true air speed is a given parameter for a segment of flight, the ground speed varies according to the wind that the aircraft is subjected to.

Let us consider an aircraft f moving at FL k from WP i on track j towards the next WP 1  i on track j  (where

j j   , 1    j j or 1    j j
) along the corresponding link (Figure 11). Let us assume that at the time t this aircraft is located at the point

with coordinates   k h , ,



. Its air speed is equal to f i v at this point (black arrow in Figure 11). The real wind speed at this point is equal to  

k h t W , ,
,    (blue arrow in Figure 11). Then, its real ground speed f V (red arrow in Figure 11) is defined by the formula:

  k j ij f i f h t W v V , , ,      , (1)
where

  k j ij h t W , ,
,    is the projection of the wind vector on the link   . To compute this time, we divide the corresponding link  

k j i , , -  k j i , , 1   .
k j i , , -  k j i , , 1   into f N sublinks and obtain 1  f N
subnodes (marked with stars in Figure 11). The ground speed

f s V at each subnode 1 ,..., 1   f N s
can be calculated via formula [START_REF] Woollings | Variability of the North Atlantic eddy-driven jet stream[END_REF]. To simplify the simulation, we suppose that this speed remains constant and equal to f s V along the link segment between subnodes s and 1  s . Let us denote the length of this link segment as s l . Then, the time

f s i t , 
necessary for aircraft f to move from subnode s to subnode 1  s is simply:

f s s f s i V l t   , . (2)
Using these time periods

f s i t , 
and given the link entry time f i t (which is the time of passing the first subnode,

f i f i t t  1 ,
), the times f s i t , at which aircraft f passes subnodes s can be successively calculated for

1 ,..., 2   f N s : f s i f s i f s i t t t 1 , 1 , ,      . ( 3
)
We thereby obtain all the necessary information to perform the simulation of aircraft flight along any link (and for the whole OTS). Note that f s i t , is the real (actual) time of arrival of aircraft f at subnode s . In addition, we consider the forecast times, noted as f s i t ,

~,

i.e. the times at which the aircraft estimates to be at subnode s using the available MF information. This time is supposed to be different from the actual time, as the exact values of real oceanic winds are not known. The

Wind Networking approach

This section describes the main idea of the WN approach. It explains how the approach was implemented in our simulations. Finally, it presents the criteria to be used for evaluating the WN benefits.

Wind Networking simulation

Let us suppose that aircraft f on route passing

the point   k h , ,  at time f t can measure the wind   k f h t W , ,
,    at this point (dark blue arrow in Figure 13, top). Measurements are assumed to be performed with some error f   . Thus, instead of the real wind function the aircraft obtains its approximation (light blue arrow in Figure 13):

    f k f k f f h t W h t W           , , , , , , ˆ. (4)
In our simulations the measurement error f 1   m/s. The value of  is set according to the statistics obtained from the previous studies [START_REF] Benjamin | Accuracy of ACARS wind and temperature observations determined by collocation[END_REF][START_REF] Cole | An assessment of the 60km Rapid Update Cycle (RUC) with near real time aircraft reports[END_REF].

Further, we suppose that aircraft can communicate the results of these measurements to other aircraft (see Figure 13,bottom). Based on the information about the winds received from preceding aircraft, the following aircraft ( g , shown in red in Figure 13) can adjust its predictions. This information (W  ˆ, red arrow in Figure 13) is expected to be more precise than MF winds (W  ~, green arrow in Figure 13), as it is more recent. Thus, the adjusted predictions should be closer to the real situation than the initial forecast. That is the main idea of WN. . To simulate this process, we divide again the link  

k j i , , -  k j i , , 1   into f N
sublinks. Using the forecast wind field W  ~, we obtain the forecast times of passing the subnodes:

f s i t , ~, 1 ,..., 1   f N s
.

Here is how we propose to improve these predictions using the WN approach.

Let us consider m aircraft

m f f f ,..., , 2 1 on
the OTS on the links neighbor to the link  

k j i , , -   k j i , , 1  
that passed the corresponding links before the aircraft f (before f s i t , ). The choice of the links, are considered to be neighbor to the current one, can vary according to user preferences. The set of neighbor links may include the current link itself, several upper/lower links and/or several northern/southern links. Here is how we define the neighbor links for most of the simulations. Let . Figure 14 displays one example of determining all neighbor aircraft (blue) of aircraft f (red) in the horizontal section of the OTS. In the vertical direction we consider only the 2 links on the adjacent FLsupper and lower, i.e.

the links   1 , ,  k j i -  1 , , 1    k j i and   1 , ,  k j i -   1 , , 1    k j i .
Aircraft f then can use the wind measurements from preceding aircraft

  k s s f s i f adj h t W , , , ˆ,    . Here f s i t , ˆ is
the aircraft predicted adjusted time of passing the subnode s , that can be calculated using formulas similar to (1)-( 3). The adjusted wind f adj W  ˆ is shown with red arrow in Figure 14 and it is supposed to be closer to the real wind (blue arrow) than the forecast wind (green arrow) at the same point.

Wind Networking evaluation criteria

The main objective of this study is to evaluate the quality of the two prediction approaches: using MF and WN. The first criterion of comparison is the error made in the time prediction. To measure these errors, we compute for each flight f and for each WP i of its route the difference between the real and the predicted times of arrival at this WP:

 f i f i f i t t e   ~
for forecast times;

 f i f i f i t t e   ˆ
for adjusted times,

and we analyze and compare the prediction errors

( f i e ~ and f i e ˆ) distributions.
Another criterion of comparison is the number of conflicts detected for the set of flights using each of the prediction approaches:

 the number of conflicts that are predicted (by each of the method) but do not happen in the reality (false alarms);

 the number of conflicts that would happen but have not been predicted (critical errors);  the errors in the conflict duration predictions (related to the time prediction errors).

The results of these comparisons are presented in the next section.

Simulation results

Several days of the years 2013-2014 were analyzed during this study. Here we present the results of simulations for 10 December 2013. The results for the other days are similar. The real and forecast wind fields used for this day are shown in Figure 6, and the valid eastbound OTS is displayed in Figure 7. Remark that the OTS passes through the area of very strong winds, which increases the uncertainty in the aircraft position determination.

Results for one track

Our first preliminary study considers a simple model consisting of several consecutive flights following one track on the same FL. For this model, several first aircraft enter every link without having preceding aircraft that already past the corresponding link, and thus, they cannot benefit from the WN approach. On the other hand, the rest of aircraft obtain complete information about the winds from all preceding aircraft for each route point. This case could very rarely happen in the reality, but it demonstrates the benefits of the WN approach under ideal conditions.

Several tests, hold with different extrapolation methods for adjusting the wind for the particular aircraft from the measurement of the preceding aircraft, revealed that using data from the latest available aircraft for the adjusted wind yields results that are as good as (and in some cases, even much better then) when applying any extrapolation method to wind values obtained from several aircraft. Thus, we focus only on these last-aircraft results in this paper.

In the example below, we consider a flight set consisting of 77 aircraft cruising on the middle OTS track (VIXUN-MALOT) at the altitude 200hPa. The flight track entry times were randomly chosen between 0000 UTC and 0500 UTC so that all aircraft exited the OTS before 0900 UTC and were separated by at least 3 minutes (according to RSS) throughout the track. First, we perform the flight simulations using the real wind field, next using the forecast wind, and finally using the wind obtained from the latest previous aircraft. For this example, the first 10 aircraft were too close one to each other to provide sufficient information (for all subnodes on links) to the following aircraft. Thus, the WN approach was only applied to the rest 67 aircraft. For each simulation, we recorded the times of passing the WPs:

f i t , f i t ~ and f i t ˆ.
Then, for the two prediction approaches, we calculated the prediction errors at each WP, i.e. f i e ~ and f i e ˆ. Figure 15 shows the two prediction-error distributions obtained. The time scale in this diagram was truncated to 30 seconds in order to focus on the mean error values.

Table 1 presents some statistics of the distributions of prediction-error absolute values. As it can be seen from this table and from Figure 15, the WN approach improves significantly the prediction of the time of passing the WPs. The average prediction error from almost half a minute (an important value for RSS) is reduced to less than 5 seconds. These encouraging results of WN application obtained for this artificial flight model give the stimulus to apply the approach to the real air traffic in NAT.

Results for real aircraft set on OTS

On the next step of our research we apply the WN approach to real air traffic data. Several data sets have been tested. In the example below we consider a set of 378 aircraft cruising between 38 airports in North America to 51 airports in Eurasia (mostly in Europe) for 4 August 2006. In order to reduce the congestion in the continental airspace, the desired entry track is assumed to be the track that is the closest to the departure airport, and similarly for the desired exit track (closest to the arrival airport) (see [START_REF] Rodionova | North Atlantic Aircraft Trajectory Optimization[END_REF] for more detail). The desired track entry time is first set to be the same as defined in the flight plan, and then it is adjusted in order to fit into the range of the OTS validity (starting from 0000 UTC). Finally, the flights are reorganized in order to avoid conflicts on links (when faster aircraft overcome slower ones), so that all remaining conflicts are conflicts on nodes.

The constructed flight set is evaluated using three different winds and the errors of time prediction of passing the WPs are compared. We have tested several extrapolation approaches to obtain the adjusted wind from neighbor preceding aircraft. Figure 16 displays the prediction error distributions for some of these approaches, i.e. when:  only the measurements received from the latest aircraft on the same link are taken into account (as in the previous subsection) (blue in Figure 16);

 the wind values obtained from the closest (in terms of distance) aircraft on the neighbor links are linearly interpolated to adjust the wind at the current aircraft position (red in Figure 16);

 the wind values obtained from the latest (in terms of time) aircraft on the neighbor links are linearly interpolated (green in Figure 16). The statistics of the distributions of the absolute values of these prediction errors are presented in Table 2. This table and Figure 16 show that using the measurements from aircraft on neighbor links improves the results; and the more recent the measurements are, the smaller the prediction errors are, and the closer to the reality the results are. Thus, in the sequel of this paper the WN approach always refers to the method based on the measurements from the most recent (latest in time) neighbor aircraft.

Figure 17 presents the comparison of the distributions of the prediction errors using MF approach (its statistics can be found in Table 2 as well) and WN approach (based on the best "latest flights" extrapolation method); the bottom diagram displays a truncated time scale in order to better demonstrate the distribution behavior near 0.

Figure 17. MF vs. WN prediction error, 378 flights

Comparison of the results from this subsection (Figure 17) with the previous results (Figure 15) shows that the WN approach seems to be less efficient in the case of several tracks and real aircraft sets: the prediction error being less than 5 seconds for one-track artificial flight set grows up to 7 seconds in the real case. Here are some explanations of this phenomenon:

 the difference in the wind evaluation for different tracks (WN is more efficient if the wind does not change much with the time);

 the time period between the consecutive flights (WN is more efficient if this period is short, which was the case for artificial flight set, but not always the case for the real one);

 the existence of previous flights (for the artificial flight set, all flights followed one another; in the real case, some flights on some links may have no preceding flights; this obliges such flights to rely on MF instead of WN, and thus, the predictions are not improved at all).

On the next stage of prediction method comparison, we concentrate on evaluating the error in conflict prediction. As mentioned above, conflicts can only happen on nodes in the considered model. First of all, we perform the flight simulation based on the three wind functions and we record the pairs of aircraft in conflict (violating longitudinal separation minimum) at each OTS node. Then, we compare the real situation with each of the prediction methods to determine the number of conflicts that were predicted correctly, the number of false alarms, and the number of errors. The first three lines of Table 3 display these numbers. The last line of Table 3 contains the absolute value of the total conflict duration time error (which is the difference between the predicted conflict duration time and the real one, considering that if the conflict is not detected its duration is 0). Figure 18 represents the distributions of conflict duration errors for both prediction approaches. The scale is truncated at 30 seconds in order to visualize better the distribution near its mean value. Table 3 and Figure 18 reveal the improvement in conflict prediction when using the WN approach. Obviously, this improvement is directly related with the accuracy of the approach in time prediction. It seems quite reasonable that the more aircraft there are the greater improvement could be achieved using WN, as more aircraft would benefit from the measurements of preceding aircraft. This hypothesis becomes even more important once remembering that the number of flights in airspace is constantly increasing. This fact encourages performing the simulations for artificially increased flight sets.

Results for increased flight set on OTS

In this subsection we present the results similar to those discussed in the previous section but for an increased flight set. In the example below we consider 1000 flights. As in any case such a flight set would be an artificial one, we chose the desired entry and exit track as well as the desired track entry time randomly. Further, these parameters were adjusted in order to eliminate the conflicts on links. We perform the simulations in the same manner as it was described in the previous subsection: we compute the prediction errors of time of passing the WPs for both approaches, the number of correctly and incorrectly predicted conflicts, the errors in conflict duration prediction.

Figure 19 presents the distribution of time prediction errors (the time scale is truncated at 1 minute). The mean absolute value of MF error in this case is 22.35 seconds; and the mean absolute value of WN error is 5.83 seconds. One can again notice the improvement in WP passing times predictions using the WN that is even greater for this case of 1000 flights in comparison to real 378 flights. Table 4 contains the number of correctly predicted conflicts, the number of false alarms and the number of critical errors, as well as the total absolute value of conflict duration prediction errors.

What is remarkable, in this case the numbers of incorrect conflict predictions as well as the conflict duration errors in average are reduced almost in 4 times when applying the WN. When observing the conflict duration prediction errors distribution (see Figure 20) one can also easily notice much greater improvement due to WN in comparison to real flight set (Figure 18). This fact reveals the absolute benefits of WN approach especially for dense traffic conditions.

Conclusion

The current paper introduces the Wind Networking approach for aircraft trajectory prediction in North Atlantic oceanic airspace. The approach is based on the possibility to use measurements from previous aircraft to adjust the forecast wind.

The paper describes the mathematical model based on which the numerical simulations were held on. During the simulations several days with different meteorological conditions were studied. The paper presents the results for one day and several flight sets that were found to be quite typical.

The results of the simulations reveal the significant decrease in the errors of prediction of times of passing the waypoints for the set of aircraft in average when applying WN approach. They also show that the number of en-route conflicts as well as their duration time are also predicted much more precisely using WN. Finally, the study performed for increased flight sets prove the special efficiency of WN in the dense traffic conditions. Thus, the WN seems to be a very promising approach to improve the flight trajectory prediction, which is especially important in oceanic airspace because of the lack of surveillance abilities. Implementing of this approach that would be possible due to new broadcast technologies will allow aircraft to know there future position in time more exactly and to determine possible conflicts more precisely in order to efficiently avoid them.

The WN brings improvement in particular in dense areas. Those areas are also the ones where a good trajectory prediction is critical. We expect in some future works to estimate the benefits of this concept in some other dense areas, like big terminal maneuvering areas.

Figure 5 .

 5 Figure 5. New surveillance/broadcast services [31]

Figure 7 .

 7 Figure 7. Eastbound OTS on 10 December 2013 We propose to model an OTS with a Nz Ny Nx  grid of WPs, where Ny is the number of OTS tracks, Nx is the number of WPs on each track, and Nz is the number of FLs for each track. The tracks are labeled from 1 to Ny starting from the most northern, the WPs on each track from 1 to Nx starting from the most western, and the FLs from 1 to Nzstarting from the lowest. Thus, the 3D position of an aircraft on track j at WP i at FL k is completely

Figure 8 .

 8 Figure 8. OTS model (horizontal section)

Figure 9 .Figure 10 .

 910 Figure 9. OTS model (vertical section)

Figure 11 .

 11 Figure 11. Simulation of flight between two WPsGiven the time

Figure 12

 12 displays the differences between real and forecast flight values. The aircraft air speed f i v is represented with black arrow, the real wind speed W  with blue arrow, while the forecast wind speed W  ~, slightly different, with green arrow. Thus, the forecast ground speed f V ~ (violet arrow) also differs from the real ground speed f V (red arrow).

Figure 12 .

 12 Figure 12. Estimation of flight between two WPs

Figure 13 .

 13 Figure 13. Adjusting wind predictions by WN Let us consider that aircraft f moves at FL k , from WP i on track j towards the next WP 1  i on track j  along the corresponding link. Being at the WP i at time f i t the aircraft needs to predict its flight until WP 1  i . To simulate this process, we

  wind, measured by this aircraft at this moment (blue arrows in Figure14).

Figure 14 .

 14 Figure 14. Adjusting wind predictions by WN Further, we can apply any appropriate extrapolation method to the corresponding measurements (5) in order to obtain the adjusted wind at each subnode 1 ,..., 1   f N s on the route

Figure 15 .

 15 Figure 15. MF vs. WN prediction error, 77 flights

Figure 16 .

 16 Figure 16. Wind extrapolation comparison

Figure 18 .

 18 Figure 18. Conflict duration error, 378 flights

Figure 19 . 4 Figure 20 .

 19420 Figure 19. MF vs. WN prediction error, 1000 flights Table 4. Conflict prediction, 1000 flights Approach MF WN Correctly predicted, Nb. 1175 1229 False alarms, Nb. 48 13 Critical errors, Nb. 70 16 Conflict duration error, min. 242.7 63.4

Table 1 . Prediction error statistics, 77 flights Approach MF WN

 1

		i e ~	f	e ˆ	i	f
	Max	215,61 26.97
	Mean	29,67 4.86
	Variance	50.83 6.04

Value (sec.)

Table 2 . Prediction error statistics, 378 flights Approach MF WN, same track WN, closest WN, latest

 2

	Mean, sec	16.36	7.48	6.87	6.40
	Var., sec	22.05	13.51	12.53	12.21

Table 3 . Conflict prediction, 378 flights

 3

	Approach	MF	WN
	Correctly predicted, Nb.	460	462
	False alarms, Nb.	11	7
	Critical errors, Nb.	12	10
	Conflict duration error, minutes 46,29 21,63

Met Office News Blog, https://metofficenews.wordpress.com/ta g/jet-stream/

Turbulence Forecast, http://www.turbulenceforecast.com/atlant ic_westbound_tracks.php

John Grimwade, Information Graphics, http://www.johngrimwa de.com/D1.html

Acknowledgements

The authors would like to thank Philippe Louyot for providing the real aircraft data files and the accompanying documents. The authors are also very grateful to Karim Legrand for productive discussions. This work has been supported by French National Research Agency (ANR) through JCJC program (project ATOMIC n°ANR-12-JS02-009-01).