Leonardo Kosloff
email: kosloff@ime.unicamp.br

Gabriele Sbaiz
email: gabriele.sbaiz@phd.units.it

Fast rotation and inviscid limits for the SQG equation with general ill-prepared initial data

Keywords: 2020 Mathematics Subject Classification: 35Q86 (primary), 35B25, 76U60, 35B40, 76M45 (secondary) SQG equation, fast rotation limit, inviscid limit, singular perturbation problem, ill-prepared data

In the present paper, we study the fast rotation and inviscid limits for the 2-D dissipative surface quasi-geostrophic equation with a dispersive forcing term AR 1 ϑ, in the domain Ω = T 1 × R. In the case when we perform the fast rotation limit (keeping the viscosity fixed), in the context of general ill-prepared initial data, we prove that the limit dynamics is described by a linear equation. On the other hand, performing the combined fast rotation and inviscid limits, we show that the initial data ϑ 0 is transported along the motion. The proof of the convergence is based on an application of the Aubin-Lions lemma.

Introduction

In this paper, we are interested in the description of the surface temperature ϑ on the ocean (see e.g. [START_REF] Cushman-Roisin | Introduction to geophysical fluid dynamics[END_REF] and [START_REF] Vallis | Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation[END_REF]). We consider the two-dimensional surface quasi-geostrophic (SQG) equation with dissipation determined by a fractional Laplacian and a dispersive forcing term in the domain Ω = T 1 × R given by:

           ∂ t ϑ + div (ϑu) + νΛϑ + AR 1 ϑ = 0 u = R ⊥ ϑ := (-R 2 ϑ, R 1 ϑ) ϑ |t=0 = ϑ 0 , (1.1)
1

where ϑ is a real-valued scalar function, u is the divergence-free velocity field and A represents the amplitude parameter for the dispersive forcing term. To define the fractional Laplacian operator Λ := √ -∆ we can adapt the Fourier transform for functions f (x 1 , x 2 ) ∈ L 1 (Ω) with (n 1 , ξ 2 ) ∈ Z × R, so that

f (n 1 , ξ 2 ) = Z×R e -i(n1x1+ξ2x2) f (x 1 , x 2 ) dx 1 dx 2 ,
and Λ is then defined by:

Λf (n 1 , ξ 2 ) = (n 2 1 + ξ 2 2 ) 1/2 f (n 1 , ξ 2 ) := | ξ| f (n 1 , ξ 2 ) .
Moreover, we have ν > 0 and R i := ∂ i Λ -1 , for i = 1, 2, are the usual Riesz transforms. The non-local dissipative term Λϑ comes from the Ekman pumping mechanism (we refer to [START_REF] Held | Surface quasi-geostrophic dynamics[END_REF], [START_REF] Lapeyre | Surface quasi-geostrophy[END_REF] and [START_REF] Pedlosky | Geophysical fluid dynamics[END_REF] for details); while the presence of an environmental horizontal gradient R 1 ϑ = ∂ 1 Λ -1 ϑ represents the advection of a large-scale buoyancy coming from the meridional variation of the Coriolis force. The parameter A > 0 is the analogue of the Rossby number which determines the typically large weight of the Coriolis term in geophysical fluid-dynamics. The mathematical analysis of (1.1) started with the work [START_REF] Kiselev | Global regularity for the critical dispersive dissipative surface quasigeostrophic equation[END_REF] by Kiselev and Nazarov, where the existence of smooth solutions in the torus was addressed.

In the case without the presence of the dispersive forcing term, the non-dissipative SQG equations (ν = 0) are the two-dimensional analogue of the 3D Euler equations in vorticity form (see for e.g. [START_REF] Constantin | Formation of strong fronts in the 2-D quasigeostrophic thermal active scalar[END_REF]), while the SQG equations are analogous to the 3D Navier-Stokes system; due to this analogy the global regularity of the SQG equations has been intensively studied in recent decades (we refer, instead, to [START_REF] Babin | Nicolaenko: Global splitting, integrability and regularity of 3D Euler and Navier-Stokes equations for uniformly rotating fluids[END_REF], [START_REF] Babin | Nicolaenko: Regularity and integrability of 3D Euler and Navier-Stokes equations for rotating fluids[END_REF] and [START_REF] Babin | Global regularity of 3D rotating Navier-Stokes equations for resonant domains[END_REF] for the pioneering studies for the Navier-Stokes equations). For an overview about these equations, we refer to [START_REF] Constantin | Long time dynamics of forced critical SQG[END_REF] and references therein. However, the presence of the dispersive forcing term makes equations (1.1) analogous to the Navier-Stokes-Coriolis (NSC) system, which is a fundamental geophysical model dealing with large-scale phenomena (we refer to [START_REF] Chemin | Mathematical geophysics. An introduction to rotating fluids and the Navier-Stokes equations[END_REF] for a more detailed discussion). The main advantage of the NSC system is that in the limit of vanishing Rossby number, the fast rotating term "produces" a stabilization effect which ensures the global well-posedness of strong solutions with large initial data, unlike the case of the Navier-Stokes equations (see also [START_REF] Cannone | Global regularity for the supercritical dissipative quasi-geostrophic equation with large dispersive forcing[END_REF] in this respect). In particular, in [START_REF] Chemin | Mathematical geophysics. An introduction to rotating fluids and the Navier-Stokes equations[END_REF] this was proved by establishing Strichartz estimates, which show how the dispersive phenomena weaken the non-linearity and stabilize NSC towards a 2D Navier-Stokes type system. Following this analogy, a similar result for the supercritical dispersive SQG equation was shown in [START_REF] Cannone | Global regularity for the supercritical dissipative quasi-geostrophic equation with large dispersive forcing[END_REF], where in (1.1) the dissipative term is represented by νΛ 2α ϑ, α < 1/2. In this context, we mention also the work [START_REF] Kosloff | Inviscid limit for SQG equation in different dispersive regimes via relative energy inequality[END_REF] (in the context of well-prepared initial data), where the main tool employed is the relative energy inequality, in order to treat the inviscid incompressible limit of the NSC system with large rotation (see [START_REF] Feireisl | Relative entropies, suitable weak solutions, and weak-strong uniqueness for the compressible Navier-Stokes system[END_REF] and [START_REF] Caggio | Inviscid incompressible limits for rotating fluids[END_REF]). Furthermore, the method used in [START_REF] Kosloff | Inviscid limit for SQG equation in different dispersive regimes via relative energy inequality[END_REF] allows to obtain the inviscid limit with fixed or no dispersion, providing an alternative to previous results [START_REF] Berselli | Vanishing viscosity limit and long-time behavior for 2D quasi-geostrophic equations[END_REF] and [START_REF] Wu | Inviscid limits and regularity estimates for the solutions of the 2-D dissipative quasigeostrophic equations[END_REF] for the unforced SQG equation.

Motivated by the previous discussion, our goal is to perform the fast rotation and inviscid limits in the more general framework of ill-prepared initial data, where strong convergence of the initial data for the limiting and target systems, or a priori structural conditions for the dispersive estimates, are not required. In this direction, we remark that, in the domain Ω = T 1 × R, the stability of the SQG system without dispersive forcing and with horizontal dissipation, was recently proved in [START_REF] Wang | Stability for the 2D anisotropic surface quasi-geostrophic equation with horizontal dissipation[END_REF]. Thus, this result shows that periodic boundary conditions allow for a decomposition into the mean flow and its oscillations, in which global existence, and indeed, stability, can be obtained in the Sobolev space H 2 (Ω), without explicit recourse to dispersive effects.

In our work, first of all, we study the regime when the rotational effects are predominant in the dynamics, keeping fixed the viscous coefficient ν > 0, i.e.

A = 1 ε , (1.2)
for a given ε ∈]0, 1]. Next, we analyse the combined fast rotation and inviscid limits where

A = 1 ε and ν = ν(ε) = ε α , α > 0. (1.
3)

The scaling for ν in (1.3) is motivated by the physical Stommel boundary layer model for the western intensification of oceanic currents, where the Ekman pumping dissipation must be taken into account. In the inviscid limit case, this scaling would typically be written with -1 < α < -2/3, but relabeling the parameter ε leads to (1.3). Note also that we disregard the frictional sublayer needed to adjust the no-slip boundary conditions (for more details, we refer to [START_REF] Desjardins | On the homogeneous model of wind-driven ocean circulation[END_REF] and [START_REF] Gérard-Varet | Remarks on Boundary Layer Expansions[END_REF]).

In order to prove our results and get the improvement to the more general ill-prepared data, on the one hand, we take advantage of the special form of test functions, belonging to the kernel of the singular perturbation operator; on the other hand, we employ the structure of the system to find compactness properties for the means of temperature ϑ ε . The special form of the test functions combined with the Aubin-Lions lemma will allow us to pass to the limit in the equations. Moreover, we need to work with periodic boundary conditions in the domain Ω, so that non-trivial test functions are still allowed by the compact support condition (see remark 3.2 after Theorem 3.1).

The strategy employed is a standard matter in the context of singular perturbations problems and it consists in the following steps:

(i) develop an existence theory, which holds for any value of ε > 0 fixed;

(ii) state uniform bounds for the family of solutions in order to extract weak limit points;

(iii) find the constraints that the limit points have to satisfy;

(iv) pass to the limit for test functions in the kernel of the singular perturbation operator.

To conclude this part, let us mention that, due to the stabilization effects of the Coriolis force, we expect convergence towards a linear equation in the case of scaling (1.2), while taking the combined scaling (1.3) we will show that the limit dynamics "collapses". This will tell us that the initial limit profile is transported along the motion.

Let us now give an overview of the paper. In Section 2 we collect our assumptions and we state our main results. In Section 3 we study the singular perturbation part of the equations, recalling the uniform bounds on our family of weak solutions and establishing constraints that the limit points have to satisfy. Section 4 is devoted to the proof of the convergence results for the fast rotation limit and the combined fast rotation and inviscid limits. Some notation and conventions. Let B ⊂ R 2 . The symbol C ∞ c (B) denotes the space of ∞-times continuously differentiable functions on R 2 and having compact support in B. The dual space D ′ (B) is the space of distributions on B. Given p ∈ [1, +∞], by L p (B) we mean the classical space of Lebesgue measurable functions g, where |g| p is integrable over the set B (with the usual modifications for the case p = +∞). We use also the notation L p T (L q ) to indicate the space L p [0, T ]; L q (B) , with T > 0. Given k ≥ 0, we denote by W k,p (B) the Sobolev space of functions which belongs to L p (B) together with all their derivatives up to order k. When p = 2, we alternately use the notation W k,2 (B) and H k (B). We denote by Ẇ k,p (B) the corresponding homogeneous Sobolev spaces, i.e. Ẇ k,p (B) = {g ∈ L 1 loc (B) :

D α g ∈ L p (B), |α| = k}. Recall that Ẇ k,p is the completion of C ∞ c (B)
with respect to the L p norm of the k-th order derivatives. For the sake of simplicity, we will omit from the notation the set B, that we will explicitly point out if needed.

We denote by T 1 the one-dimensional flat torus T 1 := [-1, 1]/ ∼, where ∼ denotes the equivalence relation which identifies -1 and 1.

In the whole paper, the symbols c and C will denote generic multiplicative constants, which may change from line to line, and which do not depend on the small parameter ε. Sometimes, we will explicitly point out the quantities that these constants depend on, by putting them inside brackets.

Let f ε 0<ε≤1 be a sequence of functions in a normed space X. If this sequence is bounded in X, we use the notation f ε ε ⊂ X.

L. Kosloff was supported by FAPESP -Brazil grant 2019/16537-0. G. Sbaiz is member of the INdAM (Italian Institute for Advanced Mathematics) group.

Moreover, the authors acknowledge Francesco Fanelli for his careful remarks that improve a lot the presentation of this work.

Setting of the SQG problem and main statements

In this section, we formulate our working hypotheses (see Subsection 2.1) and we state our main results (in Subsection 2.2).

Formulation of the problem

In this subsection, we present the rescaled SQG system with the dispersive forcing term, which we are going to consider in our study, and we formulate the main working hypotheses. The material of this part is mostly classical: unless otherwise specified, we refer to [START_REF] Kosloff | Inviscid limit for SQG equation in different dispersive regimes via relative energy inequality[END_REF] for details.

Primitive system

To begin with, let us introduce the "primitive system", i.e. the rescaled SQG system, supplemented with the scaling (1.2) presented in the introduction, where ε ∈ ]0, 1] is a small parameter. Thus, the system consists of the momentum equation and the quasi-geostrophic balance: respectively,

∂ t ϑ ε + div (ϑ ε u ε ) + νΛϑ ε + 1 ε R 1 ϑ ε = 0 (2.1)
u ε = R ⊥ ϑ ε := (-R 2 ϑ ε , R 1 ϑ ε ) . (2.
2)

The unknown is the fluid surface temperature ϑ ε = ϑ ε (t, x), with t ∈ R + and x ∈ Ω.

Remark 2.1 In the case of scaling (1.3), we replace ν with ε α , α > 0.

Initial data and finite energy weak solutions

We address the singular perturbation problem described in subsection 2.1.1, with scaling (1.2), for general ill prepared initial data, in the framework of finite energy weak solutions (see e.g. [START_REF] Cannone | Global regularity for the supercritical dissipative quasi-geostrophic equation with large dispersive forcing[END_REF]). Since we work with weak solutions based on dissipation estimates, we need to assume that the initial data satisfy the following bound:

sup ε∈ ]0,1] ϑ 0,ε L 2 (Ω) ≤ C . (2.3)
Thanks to the previous uniform estimate, up to extraction, we can argue that

ϑ 0 := lim ε→0 ϑ 0,ε , (2.4)
where we agree that the previous limit is taken in the weak topology of L 2 (Ω).

Let us specify better what we mean for finite energy weak solution (see [START_REF] Chemin | Mathematical geophysics. An introduction to rotating fluids and the Navier-Stokes equations[END_REF] for details).

Definition 2.2 We say that ϑ ε is a weak solution to the dispersive SQG system (2.1)-(2.2) in [0, T [ ×Ω (for some time T > 0) with the initial condition ϑ 0,ε , if:

(i) ϑ ε ∈ L ∞ ([0, T [ ; L 2 (Ω)) ∩ L 2 ([0, T [ ; Ḣ1/2 (Ω));
(ii) the momentum equation is satisfied in a weak sense:

for any ϕ ∈ C ∞ c ([0, T [ ×Ω), one has - T 0 Ω ϑ ε ∂ t ϕ + ϑ ε u ε • ∇ϕ + νΛ 1/2 ϑ ε Λ 1/2 ϕ - 1 ε R 1 ϑ ε ϕ dx dt = Ω θ 0,ε ϕ(0) dx ; (2.5)
(iii) the quasi-geostrophic balance is satisfied in D ′ (]0, T [ ×Ω);

(iv) the solution satisfies the following energy estimate:

ϑ ε (T ) 2 L 2 + 2ν T 0 Λ 1/2 ϑ ε (τ ) 2 L 2 dτ ≤ ϑ 0,ε 2 L 2
for all T > 0 .

(2.6)

The solution is global if the previous conditions are satisfied for all T > 0.

Under the previous assumptions (collected in subsections 2.1.1 and here above), at any fixed value of the parameter ε ∈ ]0, 1], the existence of a global in time finite energy weak solution ϑ ε to system SQG, related to the initial data ϑ 0,ε , has been proved in e.g. [START_REF] Resnick | Dynamical problems in non-linear advective partial differential equations[END_REF] (see also [START_REF] Cannone | Global regularity for the supercritical dissipative quasi-geostrophic equation with large dispersive forcing[END_REF] in this respect).

Remark 2.3

The term involving the Riesz transform R 1 does not contribute to the energy estimate (2.6), because for any s ∈ R, it holds:

Λ s R 1 θ , Λ s θ L 2 = Ω Λ s R 1 θ Λ s θ dx = - Ω in 1 | ξ| s-1 θ | ξ| s θ d ξ = -Λ s θ , Λ s R 1 θ L 2 .
(2.7)

Main results

We can now state our main results. The first statement concerns the case when the rotational and viscosity effects, with fixed ν, are predominant in the dynamics.

Theorem 2.4 For any fixed value of ε ∈ ]0, 1], assume the initial data ϑ 0,ε verifies the hypothesis in subsection 2.1.2 and let ϑ ε be a corresponding weak solution to system (2.1)-(2.2). Then, one has the following convergence property, for any T > 0:

ϑ ε * ⇀ ϑ weakly- * in L ∞ T (L 2 (Ω)) ∩ L 2 T ( Ḣ1/2 (Ω)) . (2.8)
In addition, ϑ is a weak solution to the following linear equation in R + × Ω:

∂ t ϑ + νΛ 2 ϑ = 0 , (2.9)
supplemented with the initial condition ϑ |t=0 = ϑ 0 .

The previous theorem characterize the limit dynamics of system (2.1)-( 2.2) when one consider a fast rotation regime. In contrast with [START_REF] Cannone | Global regularity for the supercritical dissipative quasi-geostrophic equation with large dispersive forcing[END_REF], in Theorem 2.4, we are able to consider data which are ill-prepared.

Instead, in the fast rotation and inviscid limits, we need more regularity on ϑ ε since we have to control Λϑ ε . In this case, the limit dynamics is trivial and tell us that the limit temperature ϑ is constant (in time). Theorem 2.5 For any fixed value of ε ∈ ]0, 1], let the initial data ϑ 0,ε ∈ H s (Ω) for s > 2. Let ϑ ε be a corresponding solution to system (2.1)-(2.2). Then, one has the following convergence property, for s > 2:

ϑ ε * ⇀ ϑ weakly- * in L ∞ T (H s (Ω)) .
(2.10)

Moreover, one deduces the relation (in the weak sense)

∂ t ϑ = 0 ,
with the initial condition ϑ |t=0 = ϑ 0 .

Remark 2. [START_REF] Bracco | Boundary layer separation in the surface quasi-geostrophic equations[END_REF] We point out that the condition s > 2 is necessary to have the embedding of Sobolev spaces H s in the space W 1,∞ of globally Lipschitz functions (see the Appendix A for more details in this respect).

Inspection of the singular perturbation

The purpose of this section is twofold. First of all, in Subsection 3.1 we recall the uniform bounds and further properties for our family of weak solutions. Then, we study the singular operator underlying to the primitive SQG equations, and determine constraints that the limit points of our family of weak solutions have to satisfy (see Subsection 3.2).

Uniform bounds

In this section we will state the uniform bounds on the sequence ϑ ε ) ε . Indeed, with the energy estimate (2.6) at hand, we can derive uniform bounds for our family of weak solutions. To begin with, we point out that, owing to the assumption (2.3), the right-hand side of (2.6) is uniformly bounded for all ε ∈ ]0, 1]. Then, one has sup

T ∈R+ ϑ ε (T ) L 2 (Ω) ≤ c (3.1) T 0 Λ 1/2 ϑ ε (τ ) 2 L 2 (Ω) dτ ≤ c for all T > 0 . (3.2)
Due to the relation (2.2) and the fact that the Riesz operator is a 0-th order operator, one gets also sup

T ∈R+ u ε (T ) L 2 + T 0 Λ 1/2 u ε (τ ) 2 L 2 dτ ≤ c , (3.3)
where the generic constant c is independent of ε > 0.

In view of the previous properties, there exist

ϑ, u ∈ L ∞ T (L 2 ) ∩ L 2 T ( Ḣ1/2
) such that (up to the extraction of a suitable subsequence) one has

ϑ ε * ⇀ ϑ and u ε * ⇀ u . (3.4)

Constraints on the limit dynamics

In this section, we establish some properties that the limit points of the family ϑ ε ε have to satisfy. These are static relations, which do not characterise the limit dynamics yet.

Proposition 3.1 Let (ϑ ε ) ε be a family of weak solutions, related to initial data (ϑ 0,ε ) ε verifying hypothesis of subsection 2.1.2. Let ϑ be a limit point of the previous sequence as identified in Subsection 3.1. Then, one deduces the relation

∂ 1 ϑ = 0 a.e. in R + × Ω . (3.5)
Proof: Let us consider the weak formulation of (2.1). We test it against εϕ where ϕ is a compactly supported test function in

C ∞ c (R + × Ω). Denoting [0, T ] × K = supp ϕ with ϕ(T, •) = 0 (and T > 0), we have -ε T 0 K ϑ ε ∂ t ϕ + ϑ ε u ε • ∇ϕ + νΛ 1/2 ϑ ε Λ 1/2 ϕ dxdt+ T 0 K R 1 ϑ ε ϕ dxdt = ε K ϑ 0,ε ϕ(0) dx. (3.6)
By uniform bounds (3.1) and (3.2), the first three integrals on the left-hand side converge to 0, when tested against any smooth compactly supported ϕ.

Analogously, thanks to the bound (2.3) on the initial data, also the term on the right-hand side of (3.6) vanishes.

Then, passing to the limit for ε → 0, we find

T 0 K R 1 ϑ ϕ dxdt = 0 , for any test function ϕ ∈ C ∞ c (R + × Ω), which in particular implies R 1 ϑ = 0 a.e. in R + × Ω .
At this point, employing the Fourier transform, one can deduce that

∂ 1 ϑ = 0 a.e. in R + × Ω . (3.7)
This completes the proof of the proposition.

Remark 3.2 Observe that the proof above relies on the fact that the compact support of test functions does not force them to be 0 everywhere, unlike the case of the domain R 2 , where the condition that the test functions are only dependent on the vertical variable implies they must vanish inside the support as well.

4 Limit dynamics

Convergence to the fast rotation limit dynamics

In this section, we will show the convergence of momentum equation (2.5) towards the linear equation (2.9) depicted in Theorem 2.4.

The uniform bounds of Subsection 3.1 are not enough for proving convergence in the weak formulation of the momentum equation: the main problem relies on identifying the weak limit of the convective term div (ϑ ε u ε ).

At this point, we rewrite the convective term in the weak formulation

- T 0 Ω ϑ ε u ε • ∇ψ dx dt = - T 0 T 1 ×R ϑ ε (u ε ) 2 ∂ 2 ψ dx dt = -2 R ϑ ε (u ε ) 2 ∂ 2 ψ dx dt , (4.1)
where the test-function ψ is defined as

ψ ∈ C ∞ c ([0, T [ ×Ω; R) such that ∂ 1 ψ = 0 . (4.2)
The previous relation (4.1) pushes our attention to find compactness properties for ϑ ε . Taking the mean in equation (2.1) and thanks to the fact that the mean of the Riesz term R 1 ϑ ε is zero, we get

∂ t ϑ ε = -div ϑ ε u ε -νΛ ϑ ε . (4.3)
To analyse the former term in the right-hand side of (4.3), we employ the paradifferential calculus (see Chapter 2 of [START_REF] Bahouri | Fourier analysis and nonlinear partial differential equations[END_REF]). Using the Bony decomposition, one can find that the term

ϑ ε u ε is in L 2 T B -1/2 2,1
and thanks to the embedding

B -1/2 2,1
֒→ H -1/2 , we argue that div ϑ ε u ε belongs to L 2 T H -3/2 . Then, recalling equation (4.3), we deduce that

(∂ t ϑ ε ) ε is uniformly bounded in L 2 T H -3/2 loc , which implies ( ϑ ε ) ε ⊂ W 1,2 T H -3/2
loc . Moreover, we already know that (

ϑ ε ) ε ⊂ L ∞ T L 2 ∩ L 2 T Ḣ1/2 .
Therefore, the Aubin-Lions lemma gives compactness of ( ϑ ε ) ε in e.g. L 2

T L 2 loc . Thus, we deduce the strong convergence (up to extraction) for ε → 0:

ϑ ε → ϑ in L 2 T L 2 loc . This implies that ϑ ε u ε → ϑ u in D ′ (R + × Ω) . (4.4)
4.1.1 Description of the limit system

With the convergence established in (4.4), we can pass to the limit in the equations.

To begin with, we take a test-function ψ (for T > 0) as in (4.2). For such a ψ, the Riesz term R 1 ϑ ε vanishes identically. Hence, the momentum equation in its weak formulation reads:

- T 0 Ω ϑ ε ∂ t ψ + ϑ ε u ε • ∇ψ + νΛ 1/2 ϑ ε Λ 1/2 ψ dx dt = Ω θ 0,ε ψ(0) dx . (4.5)
Making use of the uniform bounds of Subsection 3.1, we can pass to the limit in the ∂ t term and in the viscosity term. Moreover, our assumptions imply that ϑ 0,ε ⇀ ϑ 0 in e. g. L 2 loc . Thanks to the convergence identified in (4.4) for the convective term, letting ε → 0, we may infer that

- T 0 Ω ϑ ∂ t ψ + νΛ 1/2 ϑ Λ 1/2 ψ dx dt = Ω θ 0 ψ(0) dx , (4.6)
where we have also employed the constraint (3.5) and the special structure of the test functions (4.2). This concludes the proof of Theorem 2.4.

Proof of the convergence in the inviscid and fast rotation case

In this subsection, we consider the case when ν = ν(ε) = ε α , α > 0. In that scaling, the energy estimate reads

ϑ ε (T ) 2 L 2 + 2ε α T 0 Λ 1/2 ϑ ε (τ ) 2 L 2 dτ ≤ ϑ 0,ε 2 L 2
for all T > 0 . (4.7)

Therefore, one completely loses the uniform control on Λ 1/2 ϑ ε . In order to recover a similar bound to (3.2), we have to consider strong solutions with the following regularity properties:

ϑ ε ∈ L ∞ ([0, T [ ; H s (Ω)) and s > 2 .
The aim now is to show higher order estimates for the temperature. We have already presented in (4.7) the L 2 estimate. In order to obtain the H s control, we employ the Littlewood-Paley decomposition (see the Appendix A). Applying the dyadic blocks ∆ j to equation (2.1) (see e.g. [START_REF] Sbaiz | Fast rotation limit for the 2-D non-homogeneous incompressible Euler equations[END_REF]), we obtain

∂ t ∆ j ϑ ε + u ε • ∇∆ j ϑ ε + ε α Λ∆ j ϑ ε + 1 ε R 1 ∆ j ϑ ε = [u ε • ∇, ∆ j ]ϑ ε .
Multiplying by ∆ j ϑ ε and using the orthogonality of the Riesz term, it follows that

∂ t ∆ j ϑ ε 2 L 2 + 2ε α Λ 1/2 ∆ j ϑ ε 2 L 2 ≤ 2 [u ε • ∇, ∆ j ]ϑ ε L 2 ∆ j ϑ ε L 2 .
In particular, we get

∆ j ϑ ε L 2 ≤ ∆ j ϑ 0,ε L 2 + C t 0 [u ε • ∇, ∆ j ]ϑ ε L 2 dτ .
At this point, thanks to the commutator estimates (we refer e.g. to [START_REF] Sbaiz | Fast rotation limit for the 2-D non-homogeneous incompressible Euler equations[END_REF]) we have

[u ε • ∇, ∆ j ]ϑ ε L 2 ≤ Cc j (t) ( ϑ ε L ∞ u ε H s + ϑ ε H s u ε L ∞ ) ≤ Cc j (t) ϑ ε 2 L 2 ,
where (c j (t)) j≥-1 is a sequence in the unit ball of ℓ 2 . After summing on indices j ≥ -1, we finally derive for all t ≥ 0:

ϑ ε (t) H s ≤ ϑ 0,ε H s + C t 0 ϑ ε (τ ) 2 H s dτ . (4.8)
The scope now is finding a time T * > 0 for which the solutions are uniformly bounded (in ε) in the interval [0, T * ]. We define T * ε > 0 such that T * ε := sup t > 0 :

t 0 ϑ ε 2 H s ≤ ϑ 0,ε H s . (4.9)
Then, we deduce ϑ ε (t) H s ≤ C ϑ 0,ε H s for all times t ∈ [0, T * ε ] and for some positive constant C = C(s). Therefore, for all t ∈ [0, T * ε ] we gather

t 0 ϑ ε (τ ) 2 H s dτ ≤ Ct ϑ 0,ε 2
H s .
By using the definition (4.9) of T * ε , we finally argue that

T * ε ≥ C ϑ 0,ε H s .
In particular, this implies that there exists a time T * > 0 such that

T * := sup ε>0 T * ε > 0 .
From the definition of T * > 0 and relation (4.8), one can obtain that sup ε∈]0,1]

ϑ ε L ∞ T * (H s ) ≤ C , (4.10)
and this allows to control also the dissipative term in the relation (4.7). At this point, since those strong solutions are in particular weak solutions, one can adapt the arguments developed in subsections 3.1 and 4.1 to conclude. Indeed, letting ε → 0, one has

- T * 0 Ω ϑ ∂ t ψ dx dt = Ω θ 0 ψ(0) dx , (4.11)
for all test functions ψ defined as in (4.2). This proves the limit dynamics in Theorem 2.5.

As an immediate consequence of the first Bernstein inequality, one gets the following embedding result, which generalises Sobolev embeddings.

Proposition A. [START_REF] Babin | Global regularity of 3D rotating Navier-Stokes equations for resonant domains[END_REF] The space B s1 p1,r1 is continuously embedded in the space B s2 p2,r2 for all indices satisfying p 1 ≤ p 2 and either s 2 < s 1 -d 1/p 1 -1/p 2 , or s 2 = s 1 -d 1/p 1 -1/p 2 and r 1 ≤ r 2 .

In particular, we get the following chain of continuous embeddings:

B s p,r ֒→ W

  1,∞ , whenever the triplet (s, p, r) ∈ R × [1, +∞] 2 satisfies

	s > 1 +	d p	or	s = 1 +	d p	and r = 1 .	(A.3)

Acknowledgements

A Appendix -Littlewood-Paley theory

In this appendix, we present some tools from Littlewood-Paley theory, which we have exploited in our analysis. We refer e.g. to Chapter 2 of [START_REF] Bahouri | Fourier analysis and nonlinear partial differential equations[END_REF] for details. For simplicity of exposition, we deal with the R d case, with d ≥ 1; however, the whole construction can be adapted also to the d-dimensional torus T d , and to the "hybrid" case R d1 × T d2 .

First of all, we introduce the Littlewood-Paley decomposition. For this, we fix a smooth radial χ such that supp χ ⊂ B(0, 2), χ ≡ 1 in a neighborhood of B(0, 1) and the map r → χ(r e) is non-increasing over R + for all unitary vectors e ∈ R d . Set ϕ (ξ) = χ (ξ) -χ (2ξ) and ϕ j (ξ) := ϕ(2 -j ξ) for all j ≥ 0. The dyadic blocks (∆ j ) j∈Z are defined by 1

For any j ≥ 0 fixed, we also introduce the low frequency cut-off operator

Note that S j is a convolution operator. More precisely, after defining

for all j ∈ N and all tempered distributions u ∈ S ′ , we have that S j u = K j * u. Thus, the L 1 norm of K j is independent of j ≥ 0, hence S j maps continuously L p into itself, for any 1 ≤ p ≤ +∞.

The following property holds true: for any u ∈ S ′ , then one has the equality u = j ∆ j u in the sense of S ′ . Let us also recall the so-called Bernstein inequalities.

Lemma A.1 Let 0 < r < R. A constant C exists so that, for any non-negative integer k, any couple (p, q) in [1, +∞] 2 , with p ≤ q, and any function u ∈ L p , we have, for all λ > 0,

By use of Littlewood-Paley decomposition, we can define the class of Besov spaces. Definition A.2 Let s ∈ R and 1 ≤ p, r ≤ +∞. The non-homogeneous Besov space B s p,r is defined as the subset of tempered distributions u for which u B s p,r := 2 js ∆ j u L p j≥-1 ℓ r < +∞ .

Besov spaces are interpolation spaces between Sobolev spaces. In fact, for any k ∈ N and p ∈ [1, +∞] we have the chain of continuous embeddings B k p,1 ֒→ W k,p ֒→ B k p,∞ , which, when 1 < p < +∞, can be refined to B k p,min(p,2) ֒→ W k,p ֒→ B k p,max(p,2) . In particular, for all s ∈ R we deduce that B s 2,2 ≡ H s , with equivalence of norms:

(A.2) 1 We agree that f (D) stands for the pseudo-differential operator u → F -1 [f (ξ) u(ξ)].


Thomas Zaragoza
email: thomas.zararagozasafe@gmail.com

Yann Nicolas
email: nicolas@abes.fr

Aline Le Provost
email: le-provost@abes.fr

recognition and entity linking of contributors and their roles from statements of responsibility

Keywords: Statement of responsibility, cataloguing, access points, ner, text classification, machine learning

a classification model over a rule based approach. This pipeline is for Abes a first experience in adopting machine learning and building a generic approach with the librarian in the loop.

Abstract:

Sudoc is the french higher education union catalogue. It is run by Abes. As any large database (15 million records), Sudoc has some quality issues that can negatively impact the user experience or the database maintenance efforts, e.g. the process towards a LRM compliant catalogue.

Quality issues are diverse: data can be inaccurate, ambiguous, miscategorized, redundant, inconsistent or missing. Sometimes, they are not really missing, they are hidden, lost in some text inside the bibliographic record itself. For instance, contributor names and roles are transcribed from the document to MARC descriptive fields (statement of responsibility). Most of them have a corresponding access point that contains the normalized name and a relator code (to express the role) -optionally the identifier of an authority record. But in Sudoc, many records have contributor mentions in descriptive fields that are not identified in access points. Moreover, many access points lack a relator code. This paper will describe our efforts to extract structured information about contributors and their role from the statements of responsibility to automatically generate the following data in access points: last name, first name, relator code and optionally identifier to link to www.idref.fr, the french higher education authority file. The first step is a named entity recognition task implemented through a machine learning (ML) approach. For the recognition of names, a pre-existing generic model (from Spacy library) is employed and retrained with ad hoc data, annotated by librarians through a dedicated annotation tool (Prodigy). For roles, a model is generated from scratch. The second step is an entity linking task. The linking of contributor names is achieved with Qualinka, a logical rule based artificial intelligence framework (LE PROVOST, 2017 IFLA conference).The linking of roles is 1 Introduction Sudoc is the union catalogue of French academic and research libraries. Library catalogues are old, large and highly structured databases, created and maintained by professionals with a strong quality ethos. Recently the development of new information technology paradigms has stressed the importance of data quality: 1/ on a web of linked data, the pollution of good data by less good data is a permanent risk; 2/ predictive models and decision making algorithms require the best possible training data to optimise results and minimise the generation of additional erroneous data.

As in any large database, quality issues in Sudoc 15 million bibliographic records are diverse: data can be inaccurate, ambiguous, miscategorized, redundant, inconsistent or missing. Sometimes, they are not really missing, they are implicit, hidden or lost in some text inside the bibliographic record itself. This paper will focus on one important and interesting case.

Contributor names and roles are transcribed from the document to be recorded in MARC descriptive fields (statement of responsibility, which we will be referring to as SoR from now on). Most of them have a corresponding access point that contains the normalised name and a function relator code (to express the role) -optionally the identifier of an authority record. But in Sudoc, many records have contributor mentions in descriptive fields that are not identified in access points. Moreover, 300 000 existing person access points lack a relator code.

This paper describes our current effort to extract structured information (data) about contributors and their role from the SoRs (text). The objective is to automatically generate or correct access points containing the following data: last name, first name, function relator code and optionally identifier to link to www.idref.fr, the french higher education authority file.

Our effort can be decomposed into two main parts: 1/ creating and evaluating a Name Entity Recognition (NER) model to extract Person and Role entities from the SoRs; 2/ linking the extracted role to the UNIMARC controlled vocabulary of roles via text classification. As this project is still a work in progress, the last part will give an overview of further work to be done.

Incomplete bibliographic records: Missing access points.

Given a UNIMARC bibliographic record with the following statement of responsibility (SoR):

200 1 $aHawking$fStephen Finnigan, réalisation$gStephen Hawking, Stephen Finnigan, Ben Bowie, scénario$gJoe Lovell ; Tina Lovell ; Arthur Pelling [et al.] acteurs B200$f and B200$g UNIMARC fields encode the SoR as found on the document (mainly the title frame as a source for a motion picture or the title page for a book etc.). The SoR transcribes, "records", the original text found on the document, with minimal structuration: a separate SoR per role/function. Hence, a unique SoR can mention more than one person: $gStephen Hawking, Stephen Finnigan, Ben Bowie, scénario

The following UNIMARC fields deal with access, not description. This highly structured data comes from the intellectual analysis of the document by the cataloger, not the transcription of the title page: Each line is called an access point, and refers to a unique person.

In 701 1$3241177421$aLovell$bTina$4005, the subfields respectively refer to the linked authority record 1 , the last name, the first name and the UNIMARC role code 2 of the person relative to this document.

It is very easy for a person reading the previous SoRs to count the number of different persons mentioned (6) then to compare to the amount of access points (5) and come to the conclusion that this record is lacking access points. But this conclusion is not trivial to reach automatically.

In the process of encoding the 5 access points, the cataloguer extracted person names from the 3 SoR fields, but seems to have dismissed or forgotten one person (Ben Bowie). We want to rely on machines to extract the missing person, and to predict which precise role (function) Ben Bowie has. 300 000 author access points lack a relator codelacunae that we will try to fill. Missing access points in the Sudoc database are, by construction, more difficult to find: this task would require to know how many different names are mentioned in the SoR. That is precisely one of the objectives of the project, which we can accomplish using a named entity recognition model. 1 https://abes.fr/es/reseaux-idref-orcid/le-reseau/ 2 https://www.ifla.org/files/assets/uca/unimarc_updates/BIBLIOGRAPHIC/u_b_appb_update2 020_online_final.pdf

'Person' and 'Role' entity extraction (NER)

Named Entity Recognition (NER) is an application of Natural Language Processing (NLP) on large amounts of unstructured text to extract entities. We've decided to apply open-source models offered by spaCy3 just for this purpose.

The model we use is 'fr_core_news_lg', that has been trained on 7200 high quality, hand annotated french articles on wikipedia (WIKINER) and over 3000 french sentences (French Sequoia). It can be used to extract persons, organisations and locations but not roles. For the latter, a new model is created and trained from scratch.

By applying the initial model on the statement of responsibility:

Stephen Finnigan, réalisateur we obtain:

And whilst it may offer promising results in a broad context, this result doesn't translate as much in a bibliographic context. Statements of responsibility are not written as naturally as the articles and sentences this model was trained on. Moreover, this model ignores the Role entity. So it was important to be able to retrain this model in a bibliographical context, and for this we had to annotate a large quantity of bibliographic records.

To do this, we have used Explosion's4 annotation tool: Prodigy. Prodigy is an annotation tool that allows us to annotate a greater amount of data faster and easier, with spaCy's models pre-built-in for training and retraining models. While we can retain spaCy's capacity for recognizing people, it has no capacity to recognize roles.

After annotating 10 000 bibliographical records using Prodigy, we were able to retrain spaCy's model to better extract Person entities and begin extracting Role entities.

Evaluation

A new model can be evaluated using three metrics: the precision, the recall and the accuracy. These metrics can be visualised using a confusion matrix.

Figure: Confusion matrix example

The precision is the chance for a prediction to be true. It is calculated by dividing the amount of true positives by the sum of true positives and false positives. The recall is the chance for an entity of a certain class to be predicted as such. It is calculated by dividing the amount of true positives by the sum of true positives and false negatives. Finally, the overall accuracy of a model is the sum of true positives and true negatives divided by the sum of all four.

To return to our previous problem, after retraining our model, we have obtained the following results: tackle our next objective while still being able to further retrain our model, given additional annotated data.

To reiterate, we can now use a NER model to analyse statements of responsibility to detect the people mentioned and the keywords that represent their role in the creation of the catalogued document.

4 Linking role keywords to their role code using a text classification model

Now we want to link role keywords to their controlled role code. To illustrate the issue, let's once more rely on the following example:

Stephen Finnigan, réalisateur

There are around 155 different relator codes. For a human it is quite easy to link 'director' ('réalisateur') in the context of a movie to the function relator code 300 -Movie director. But there can be difficulties clearing ambiguity without enough context. A director could also be an 632 -Artistic director, 727 -Thesis Director etc.

To automatically remove this ambiguity, we have first approached an Entity-Linking model.

An entity-linking model consists in linking entities of a text to a controlled ontology. The problem we encountered with this approach was that it relied too much on a rules-based system, which could produce excellent results, but also required a larger quantity of work that was not viable within the context of this project.

So instead we have opted for a text classification approach. Text classification is a machine learning technique that assigns a set of predefined categories to open-ended text through various variables like word count frequencies and other predefined features. This would require no set of rules. The human work is limited to choosing features, extracting data, creating training sets, and choosing an algorithm that provides the best result for our problem.

Initially, we have decided to provide the following features: the role keywords 6 , the document content type 7 , the position of the mention of responsibility 8 and the document type 9 (If thesis or not).

Creation of the training dataset

To get the previous features and create a training set, we extract from the Sudoc database the entries for the UNIMARC fields 200$ (SoR), 70X$ (access points), 181$c (content type), 608$3 (ID of the authority record of the form/genre of the document), 105$a (textual resources types) and 503$a (Form title).

With our NER model we have extracted from the SoR the role keywords and the person names. From which we will pair these names with the persons in the access points to provide a relator code 'answer' for the training of our model. Here are the main steps of this process:

Once the raw data is extracted we note the position of each SoR, deduce whether the document is a thesis or not. We obtain the following entries: The most complex step is executing our NER model on the SoRs to extract persons and roles keywords Our approach has limits. Since it is vital that the training set is as accurate as possible we must ignore any ambiguity encountered that we are unable to remove. For example, if this bibliographic record was correctly filled out, Stephen Finnigan would have the following roles: 300 and 690. And in such case, it would be impossible for us to be able to distinguish between the associated keywords and function relator codes in the following entry: We first implemented the KNN model as an initial test to study the feasibility of our approach before branching out and trying various different algorithms and choosing the best suited, because the KNN algorithm is fairly simple to comprehend and makes for a nice introduction to the world of machine learning and decision making algorithms.

The entries of our training set are encoded into numerical values that can be computed by our model. Each training entry is then a point in a universe. For example, let's say we train a KNN model with k=3 and two classes. Author is in red and Illustrator in blue.

Figure: Prediction by KNN model with k=3

When attempting to predict the classification of a new point (grey in this case), we will look for its 3 nearest neighbours and predict the majority. In this case, the role Illustrator in blue.

Evaluation

After training our model, we can visualise the results in the following The greatest discrepancy is in the 070 -Author column of predictions. This class has a bad precision of 0.444 which is due to the nature of said class. Being the most common and principal role in the Sudoc, often in first position and not introduced by a role keyword, we implicitly arrive at the conclusion that said person is the author. Thus this lack of role keywords make this class a sort of 'collective bin' for all other entries that either lack role keywords (due to either a mistake by the NER model or initial cataloguing) or are simply too different to the average entry of its class. In practice since the class 070 -Author is a majority role (see pie chart on distributions of roles per access point), the fact that class 070 -Author is the de facto prediction for sparse entries is the least damaging to overall efficiency. And the decent recall and precision scores for other classes assure that the more specific and least common classes are correctly predicted. As mentioned previously though, there are over 100 roles that our model needs to be able to differentiate and predict between. Considerably more than 13 roles. Thus it is unreasonable to believe our current approach is not flawed. Indeed if we train another model, this time differentiating between 35 classes, we obtain the following result: We have outlined the cells with a score inferior to 0.7 and we can deduce this time that there are many more 'collective bins'. This is due, in part, to the similarities between these roles that did not exist in the previous. If we chart out (simplified) the mispredictions: We notice the apparition of two specialised 'collective bins', the classes 350 -Engraver and 673 -Director of the research team. In this case, these classes represent respectively the artistic and academic classes. The existence of these collective bins shows there is not enough features to distinct between similar roles. We have to diagnose the mispredictions and determine if and how we can better our results.

-Further work

Firstly we wish to optimise and improve on the results achieved and work accomplished.

There are over a hundred function relator codes and we were able to implement a text classification model able to predict only for 33 function relator codes. And we were able to obtain good enough results to justify an eventual implementation of our models on a grand scale for only 13 function relator codes among the 33 function relator codes.

To achieve this, we plan to vastly increase the amount of data we can allocate to our text classification model's training. As we are currently constrained by the issue of keeping our classes balanced, we can both increase the amount of bibliographical records extracted and used to automatically create training entries and/or manually extract records including the lacking function relator codes to create specific training entries.

We also plan to increase the quantity of features that a model may rely upon when making its decision to increase the quality of the prediction and better differentiate between similar function relator code. We are currently working hand in hand with a team of cataloguing experts to analyse the false predictions and explore new potential features.

Another possible solution would be to create different models for different needs. To elaborate, we could create and train multiple models that are more specialised to a certain category of documents rather than using a large all-encompassing model. This would allow us to both reduce class constraints when it comes to creating balanced training sets but also allow our models to be able to better differentiate between similar role codes as it will be able to focus on the finer details and differences.

We can also implement different classification algorithms to find one better suited to our needs and not limit ourselves to the KNN algorithm.

And we may also implement an additional decision making algorithm above the current prediction model to increase the accuracy. The aim would be to analyse the score given by the initial prediction model and only act upon it at a certain threshold. To present this idea in simpler terms, we would only take into consideration the prediction if it is 'sure enough' and ignore it if it is 'unsure'. This would increase precision at the cost of recall, but this is a sacrifice we would be willing to take given the nature of the project. This would also give us the option to accept multiple predictions instead of the initially single one offered.

We would also want to take into consideration the hierarchy between roles. For example, an engraver ( 350) is a more specific type of author ( 070).

Finally we have to implement the generation of missing access points. It requires to differentiate between first and last names to 1/ Use the logical rules based system tool Qualinka10 , to link the person names to their idRef identifier. 2/ If it cannot be linked, we will still need to differentiate first and last names when creating the UNIMARC zones and subzones for the access point.

-Conclusion

The objective of this project is to teach models to analyse the text of statement of responsibility to automatically generate missing access points for contributors -or adding the function relator code when missing in an existing access point.

It first requires the extraction of two types of entities: persons and roles. The standard Spacy model used for this NER task gives excellent results (precision > 0.99 and recall > 0.92) after being re-trained on 10 000 manually annotated records. This annotation task, critical but often time consuming, was efficiently achieved with the help of a dedicated annotation web tool. Within this comfortable environment, AI assisted librarians can produce more reliable training data which will result with further improvement of the AI.

When the persons and roles are extracted and paired, the next arduous task is to find the right function relator code among UNIMARC relators. Our first results with KNN as classification algorithm will be completed and hopefully enhanced in various ways: more data, more features, more models (specialised), more algorithms, … and more librarians to analyse the predictions and work with the data scientist.

This project is still a work in progress. If the end result is satisfactory, it could be used in production to create new access points or check if the existing ones are correct. But no matter the final conclusion of the project, it will have taught Abes a lot: 1/ As a bibliographic agency and a massive "data steward", Abes has the duty to quickly adopt the machine learning approaches to fulfil its traditional missions ; 2/ with the help of efficient tools to annotate and prepare the data and with librarians in the loop, we can achieve real progress both in terms of data quality and human resource development, two main issues for libraries.

  3028590295 $aHawking $bStephen $4690 701 1 $3241177286 $aLovell $bJoe $4005 701 1 $3241177421 $aLovell $bTina $4005 701 1 $3241177588 $aPelling $bArthur $4005

Figure :

 : Figure: Confusion matrix of our spaCy model retrained on 10 000 records

  . The keywords are linked to the appropriate persons using a simple Person* + Role* / Role* + Person* pattern: Afterwards we make a string comparison of the Person entity and the persons in the access points taking in consideration abbreviations, composite names and eventual typos using the Levenshtein distances between the two names. Given the NER model has correctly extracted the entities, the extracted name and access point name have been paired, we can create the following training set:

Figure :

 : Figure: Pie chart of the distribution of function relator codes per access point for all records

Figure :

 : Figure: Prediction relation graph between the function relator codes

  Directeur' could be associated with the role code 300 or 690 and it requires a decision to be made to choose between two… which is the objective of our model. A conundrum as we can not use the model's predictions to create its own training set. At least not yet.The inability to differentiate between multiple role codes is quite limiting and a root for discrimination for the training of our model as some roles have a higher chance of appearing in conjunction with others than by itself. So these roles tend to have less training entries.

	Bibliographic	Keywords	Position	Content type	Thesis	Label
	Record ID
	236018256	Directeur,	0,1	tdi	False	300,69
		scénario
	Nonetheless we were able to create a training set of 1000 entries per function for 13 of the
	most common functions: 005 -Actor , 065 -Auctioneer, 070 -Author, 100 -Original
	Author, 230 -Composer, 300 -Movie director, 340 -Scientific publisher, 365 -Expert, 440
	-Illustrator, 651 -Publication director, 727 -Thesis director and 730 -Translator. And
	created then trained a KNN algorithm model

'

We are limiting ourselves to 1000 entries as whilst we may have had hundreds of thousands of annotated training entries of certain classes (like the role code 070 -Author) we only have thousands of entries of minority classes (like the role code 300 -Composer). To avoid creating a biassed predictive model we have to create a balanced training set. The easiest and safest method in our case was 'downscaling', which entailed creating a training set with an equal amount of entries per role code. That amount is the highest possible while still remaining inferior to the amount of training entries from the lowest count class.

K-Nearest Neighbour (KNN)

"If it looks like a duck, if it sounds like a duck, chances are… it's a duck."

  table and confusion matrix:That shows an overall decent prediction model, given the limited amount of data used to train and test it.

	Label	Precision	Recall	F1-Score	Support
	005	0,99	0,97	0,98	224
	065	0,93	0,78	0,85	211
	070	0,44	0,91	0,6	196
	100	0,96	0,81	0,88	217
	230	0,99	0,98	0,98	179
	300	0,92	0,98	0,95	191
	340	0,95	0,61	0,74	201
	365	0,97	0,9	0,93	181
	440	0,96	0,78	0,86	197
	651	0,81	0,6	0,69	58
	690	0,94	0,96	0,95	206
	727	0,95	0,99	0,97	207
	730	0,99	0,89	0,94	205
		Figure: Confusion Matrix for KNN model, K=6 , 13k training

https://spacy.io

Team behind development of spaCy.

Le Provost, Aline and Nicolas, Yann. "IdRef, Paprika and Qualinka. A toolbox for authority data quality and interoperability"ABI Technik, vol. 40, no. 2, 2020, pp. 158-168. https://doi.org/10.1515/abitech-2020- 2006

Acknowledgments

The authors are grateful to Pascal Poncelet ( 003), Abes colleagues who helped us to extract and analyse data, and Stephen Finnigan.


Zhengyi Wang
email: zhengyi.wang@enac.fr

Daniel Delahaye
email: daniel.delahaye@enac.fr

Jean-Loup Farges
email: jean-loup.farges@onera.fr

Sameer Alam
email: sameeralam@ntu.edu.sg

Route network design in low-altitude airspace for future urban air mobility operations A case study of urban airspace of Singapore

Keywords: Urban Air Mobility, Route Network Design, Urban Airspace, Air Traffic Management

The rapid growth of the metropolis leads to the proposal of alternative solutions, including the concept of Urban Air Mobility (UAM). Automated and highly-integrated UAM networks have been proved to have great advantages to handle UAM traffic flow with intense density and complexity in the near future. This research addresses designing UAM route networks in low-altitude airspace to minimize noise impact and maximize the efficiency and safety of UAM operations. Singapore's urban airspace is selected for the case study. On the basis of the open-source data of Singapore, the UAM network is designed as a grid-based multi-layer route network that supports twoway traffic. The topology features of the route network are analyzed. To provide alternative travel options for UAM traffic flow, we search for feasible routes between Origin-Destination (OD) pairs by solving k-Shortest Path with Diversity (KSPD) problem that minimizes link costs in terms of noise impact, safety and efficiency. The resulting feasible routes can potentially reduce airspace complexity and can be used for air traffic assignment. In addition, the impact of different parameter settings for link costs on UAM services is explored.

INTRODUCTION

Urban Air Mobility (UAM) is an air transportation concept that supports passenger or cargo-carrying air transportation services in and around urban environments [START_REF]Urban Air Mobility (UAM) concepts of operations v1.0[END_REF]. European Union Aviation Safety Agency (EASA) claimed that by 2024-25, UAM may be a lived reality in Europe [START_REF]Study on the societal acceptance of urban air mobility in europe[END_REF]. To categorize the expected evolution stages of the UAM transportation system, NASA has developed a framework for UAM Maturity Levels (UML) [START_REF] Goodrich | Description of the nasa urban air mobility maturity level (uml) scale[END_REF]. During the intermediate and mature state of UML, the traffic density of UAM vehicles becomes intense. For high-density and complexity UAM operations in densely populated metropolitan areas, a centralized UAM transportation system with highly-integrated and automated route networks has great advantages [START_REF]Urban Air Mobility (UAM) concepts of operations v1.0[END_REF][START_REF] Goodrich | Exploratory analysis of the airspace throughput and sensitivities of an urban air mobility system[END_REF].

Significant progress has been made to design UAM route networks in urban airspace. FAA suggested adopting existing helicopter routes to construct route networks for early-stage UAM operations [START_REF]Urban Air Mobility (UAM) concepts of operations v1.0[END_REF] . Tang et al. generated a 3D route network for all pairs of vertiports by using the visibility graph to avoid obstructions [START_REF] Tang | Automated flight planning of high-density urban air mobility[END_REF]. Some studies focus on UAM hub-and-spoke network design [START_REF] Wu | Integrated network design and demand forecast for on-demand urban air mobility[END_REF][START_REF] Willey | A method for urban air mobility network design using hub location and subgraph isomorphism[END_REF]. Hong et al. designed a delivery network to support commercial stand-alone drone delivery service in an area with obstacles [START_REF] Hong | A range-restricted recharging station coverage model for drone delivery service planning[END_REF].

Many studies have only focused on on-demand transportation services. Since on-demand services are often represented by a fully-connected graph [START_REF] Willey | A method for urban air mobility network design using hub location and subgraph isomorphism[END_REF], the network size has to be limited in some studies. In addition, very few studies consider environmental impact in network design. The lack of integration with real-world data is also a shortcoming of some studies.

To fill the aforementioned research gap in UAM route network design and to minimize the impact on existing air traffic management systems, this paper proposes a methodology to design a UAM route network in low-altitude urban airspace in the presence of obstacles and hazardous airspace. Using a variety of open-source data, this approach is applied as a case study for parcel delivery services using Unmanned Aerial Vehicles (UAVs) as UAM vehicles in Singapore's urban airspace. Firstly, The UAM route network is designed as a grid-based network that avoids obstacles and airspace within which the aircraft is prohibited. Then, the link cost is defined in terms of noise impact on populations, airspace safety, and flight efficiency. Unlike most researches that only calculate the shortest path, we formulate a k-shortest path problem with diversity to select feasible routes between Origin-Destination (OD) pairs that minimize the route cost. The feasible routes with low similarity will provide more travel options, which can be used as pre-computed routes for traffic assignment models [START_REF] Dafermos | The traffic assignment problem for a general network[END_REF][START_REF] Wang | Air traffic assignment for intensive urban air mobility operations[END_REF] to support high-density and complexity flowbased operations.

DATA SOURCE AND DATA DESCRIPTION

To model the static obstacles in the urban airspace of Singapore, we use ALOS World 3D-30m (AW3D30) data [START_REF]Alos global digital surface model[END_REF]. It is a global Digital Surface Model (DSM) that represents land terrain (natural and built features) with approximately a resolution of 30 meters. The ALOS World 3D data that covers Singapore is visualized in Figure 1. To identify airspace that is potentially hazardous to UAM operations, the airspace classification information of Singapore [START_REF]Airspace classification information of singapore[END_REF] is used in this study. These airspace types are promulgated as follows [START_REF]Air Navigation (Protected Areas) (No. 2) Order[END_REF][START_REF]Air navigation regulations (101 -unmanned aircraft operations)[END_REF]: In addition, the location of vertiports is collected from a study that determined the optimal distribution of vertiports for Singapore considering both costs and society value [START_REF] Zeng | Future demand and optimum distribution of droneports[END_REF]. Some locations of vertiports are modified to avoid hazardous airspace and terrain. Some delivery sites are randomly generated in the usable airspace. The location of vertiports and delivery sites are depicted in Figure 2.

•
To measure the noise impact of the UAM network on populations in different regions, the Singapore residents (as of June 2021) by planning area [START_REF]Singapore residents by planning area / subzone, age group, sex and type of dwelling[END_REF] and the boundary of planning area [START_REF]Master plan subzone boundary (no sea)[END_REF] data are collected. Singapore is divided into 332 subzones, and the total resident population in each subzone is illustrated in Figure 3. Note that in some subzones, the population may be nil or negligible.

METHODOLOGY

In the low-altitude urban airspace of a metropolitan area, the problem addressed here aims to design a multi-layer route network connecting Origin-Destination (OD) of UAM operations that optimizes the airspace usage, flight efficiency, and noise impact on residents, while keeping away from terrain and hazardous airspace. Based on the aforementioned datasets, we firstly determine the minimum spatial area to design the network and partition the area into grids. The grid points and air routes that have no interaction with the terrain and hazardous airspace are considered nodes and links, respectively. We search for the k-minimum cost path between each OD, where the link costs take into account airspace usage, noise impact, and flight efficiency.

A. Urban airspace usage and obstacles

According to the literature [START_REF] Cho | How to assess the capacity of urban airspace: A topological approach using keep-in and keep-out geofence[END_REF], several terms that describe the availability of airspace are firstly defined. The raw availability of urban airspace is represented by free airspace A f , which is the airspace free of static obstacles Ω. Free airspace can be further divided into unusable airspace A u and usable airspace A s . Unusable airspace includes prohibited areas, protected areas, and airspace affected by geofence associated with obstacles. The rest of the free airspace is denoted as usable airspace. To avoid interfering with other air traffic, we create network layers on altitudes 150ft, 175ft, 200ft, 225ft, and 250ft. In different layers, we consider the terrain above the current altitude as static obstacles. To maintain safety separation between UAVs and obstacles, we apply geofence on terrain to impose geographical boundaries around the perimeter of obstacles as a virtual barrier. We define the minimum keepout distance as d k =30m in accordance with the specification of most countries, as illustrated in Figure 4.

B. Grid-based network

Firstly, World Geodetic System (WGS) is converted to Projected Coordinate System (PCS) EPSG:3414 for Singapore. To reduce the search space and the network construction cost, we calculate the Axis-aligned Minimum-Area Bounding Box (AMABB) of vertiports and delivery sites as the minimum spatial area to design route networks and discretize the space with waypoints. We use the method of rotating calipers [START_REF] Toussaint | Solving geometric problems with the rotating calipers[END_REF] to calculate AMABB. The idea is to calculate the convex hull, then rotate by each edge orientation and compute the bounding rectangle. The minimum intersection area of the bounding rectangle and land of Singapore is the minimum area for delivery. Then, the area is discretized with grids starting from a corner point. To reduce the computation time while providing the right amount of node and link options, the size of discretization is set as d g = 2000m. The diagram of these steps is illustrated in Figure 5.

The obstacles with geofence in the delivery area at different altitudes are shown in Figure 6. The urban airspace especially below 40m contains a large number of geographical constraints, which seriously limits the connections between nodes in the route network. This in turn causes bottlenecks of congestion, especially at intersections. which reduces overall airspace capacity and performance. As a remedy, the route network has a high level of flexibility and connectivity in the airspace at higher altitudes, where few or no static obstacles still exist.

At each altitude, the grid points that do not collide with obstacles, the points that have the same longitude and latitude as vertiports and delivery sites are selected as potential waypoints. Any pair of potential waypoints with a distance not greater than √ 2d g and not colliding with unusable airspace A u and obstacles Ω are connected as a link. In addition, the vertiports and delivery sites, and all waypoints directly above them are vertically connected as links. All these links form the set of link E. All nodes that connect links in E form the set of node V. The resulting route network is represented as G = (E, V).

C. Constraints and cost for feasible paths

To support high-density and complexity UAM traffic flow in route network G, the candidate paths P = {P w |w ∈ W} should be determined for a given set of OD pair W. We search for candidate paths that maximize flight safety and efficiency, minimize the noise impact on populations in the nearby subzones, and satisfy operational constraints in G.

The length of path d p should not exceed the range of UAV, denoted as R max :

d p < R max , p ∈ P (1)
The flight time should be limited to the endurance of UAV:

e∈p d e s e < T max , p ∈ P (2)
where s e is the average speed of UAM traffic flow on link e. Note that, the maintenance and battery charging (or refueling) is assumed to be available both at vertiports and delivery sites.

If not, the range and endurance of the UAV should be cut by half to consider the return trip.

It is noteworthy that some part of usable airspace may still be at risk, including aerodromes A a , danger areas A d , and prohibited areas A p . The link traversing these areas is characterized by δ e,n and penalized by a coefficient α n , n ∈ {a, d, p}. If link e traverse hazardous airspace n, δ e,n = 1, otherwise δ e,n = 0. According to the hazard level of airspace [START_REF]Air Navigation (Protected Areas) (No. 2) Order[END_REF][START_REF]Air navigation regulations (101 -unmanned aircraft operations)[END_REF], α a , α d , and α p is set to be 0.2, 1, and 1, respectively. Then, we define the cost for crossing hazardous airspace C o,e for link e ∈ E based on the airspace type:

C h,e = (α a δ e,a + α d δ e,d + α p δ e,p ) (3)
The cost of flight efficiency is calculated based on link length d e and the type of link:

C f,e = (ϕ h δ e,h + ϕ u δ e,u + ϕ d δ e,d )d e (4)
where d e is the length of link e, δ e,h = 1 if link e is horizontal, otherwise 0. δ e,u and δ e,d similarly. ϕ h , ϕ u , ϕ d are efficiency coefficients for horizontal, upward and downward links, respectively. The coefficient is much larger for vertical links than horizontal links because crossing altitudes will introduce traffic complexity in another dimension. The link for climbing is more penalized than the link for descent as it is more energy-consuming. Considering the distance gap between horizontal and vertical links, to make their costs of flight efficiency comparable, we set ϕ h = 0.005, ϕ u = 1 and ϕ d = 0.5. Another focus is to mitigate the noise impact of the designed link to populations in each subzone of Singapore. We use Sound Pressure level (SPL) to measure the pressure level of noise produced by UAVs [START_REF] Scozzaro | Noise abatement trajectories for a uav delivery fleet[END_REF]. We assume the UAV of height h produces the same reference noise of L h at the ground point directly below, as shown in Figure 7. In addition, we denote the estimated UAV density on link e is N e . Based on the inverse proportional law, the SPL L r at the distance r from the point directly under UAV of density N e can be simplified as:

L r = L h + 10 log 10 N e h 2 h 2 + r 2 (5)
The noise impact of the UAV traffic flow on a link e to a nearby subzone S can be represented as the integral of the SPL from a point distributed on the link to a point distributed in the subzone:

Le = e f e (b) S f S (a)L dist(a,b) dadb (6)
where dist(a, b) is the euclidean distance between point a and b, f S (a) is the Probability Density Function (PDF) of a point a in polygon S, f e (b) is the PDF of a point b in line segment e. The details for calculating Le is given in Appendix A.

We search for subzones with polygonal shapes in the vicinity of each link modeled as a line segment. The total noise impact from a link e to all nearby subzones is defined as:

C n,e = dist(S,e)≤dt

P S Le (7)
where P S is the population in subzone S, d t is the distance threshold. To avoid significant effects by other impact factors such as wind, d t is set as a short distance 2000m. dist(S, e) is the distance from subzone S to link e and is defined as:

dist(S, e) = e f e (b) S f S (a)dist(a, b)dadb (8)
The calculation of equation ( 8) is the same as equation ( 6). The total cost for a link e is given as:

C e = w h Ch,e + w f Cf,e + w n Cn,e , e ∈ E (9)
where w h , w f , w n are weighting parameters between 0 and 1 and satisfy w h + w f + w n = 1, Cn,e , Ch,e , Cf,e are normalized costs. The prior scaling of the objectives are mandatory to make them comparable. Min-max normalization is used in this study:

C = C -min(C) max(C) -min(C) (10)
The path cost is the sum of the costs of links that belong to this path:

C p = e∈p C e (11)

D. Finding k-shortest path

For a route network G in a metropolitan area, the number of paths connecting each OD pair may be very large. To reduce the computational burden and potentially decrease the airspace complexity, a natural idea is to reduce the number of candidate paths between each OD pair by searching the k-Shortest Paths (KSP) based on the path cost defined previously. Denote the acyclic paths satisfying constraints (1-2) as P, then the problem is to find for each OD w the set of k paths P w satisfying ∀p ∈ P/P w , ∀q ∈ P w , C p ≥ C q .

However, in KSP problems, it is very likely that some of the shortest paths are highly similar. They share a large part of subpaths, and the resulting paths provide low flexibility and few alternatives for air traffic assignment. To this end, we search k-Shortest Paths with Diversity (KSPD). Since KSPD has been proved to be NP-hard [START_REF] Liu | Finding top-k shortest paths with diversity[END_REF], we use a greedy algorithm to approximate the KSPD problem. A similarity metric is defined as follows:

sim(p, q) = C Ep∩Eq 2C p + C Ep∪Eq 2C q (12)
where E p is the set of links on path p.

If the similarity between two paths is higher than a threshold, they are considered similar paths and cannot exist at the same time. The pseudocode of this algorithm is given in Algorithm 1. The idea is to insert the next shortest path that is sufficiently dissimilar to all shortest paths in the result set. This algorithm has been shown to have a high approximation ratio to the exact solution in real-world applications [START_REF] Liu | Finding top-k shortest paths with diversity[END_REF]. The similarity threshold τ p is chosen to be 0.6.

RESULT AND DISCUSSION

A. Experiment setting

According to acoustical measurement data of UAV in [START_REF] Schäffer | Drone noise emission characteristics and noise effects on humans-a systematic review[END_REF], we roughly set L 150ft = 80dB, L 175ft = 78dB, L 200ft = 76dB, L 225ft = 74dB, and L 250ft = 72dB. Referring to [START_REF] Troudi | Sizing of the drone delivery fleet considering energy autonomy[END_REF], the range and endurance of UAV are set as R max = 25km and T max = 1h. The average speed of UAM traffic flow on horizontal links and vertical links are set to be 15m/s and 2m/s. For each OD pairs, we select the number of shortest paths k = 20, and the estimated density of UAV on all links is the same. The framework is implemented in Python 3.7 on a laptop with Intel i7-8750H CPU and 32GB DDR4 RAM. for w ∈ W do To compare the horizontal network in different altitude layers, by removing vertical connections, some basic topology features are calculated on the basis of network theory. These topology features include the number of nodes, number of links, average degree, average clustering coefficient [START_REF] Saramäki | Generalizations of the clustering coefficient to weighted complex networks[END_REF], and degree of assortativity [START_REF] Newman | Mixing patterns in networks[END_REF]. Their values for different altitude layers are summarized in Table I.

As is shown in Table I, the network has more number of nodes and links at higher altitudes, because higher altitude layers have fewer obstacles. The result of the average degree indicates that each node has more connections to neighbors as altitude increases. The growth of the average clustering coefficient with altitude signifies that nodes tend to cluster together as altitude increases. Except for altitude layer 150ft which is non-assortative, the degree of assortativity for other network layers is positive, which means that high degree nodes are more likely to attach to high degree nodes, that is, these networks are assortative.

C. Link costs

The normalized link cost for crossing hazardous airspace Ch,e , of flight efficiency Cf,e , and of noise impact Cn,e are calculated and represented in Figure 9, 10, and 11, respectively.

According to Figure 2, since hazardous airspace occupies a large amount of Singapore's airspace, most of the links in the route network have to traverse hazardous airspace. This fact is consistent with Figure 9. Nevertheless, airspace within 5km of aerodromes takes up the largest proportion of hazardous airspace, so most links have a low cost for traversing hazardous airspace. There are some areas where different types of hazardous airspace overlap, where link costs are significant. In terms of link cost of flight efficiency shown in 10, for horizontal links, longer links are less efficient and have higher costs. Though vertical links are heavily penalized, it can be seen that the costs of flight efficiency for vertical links are still less than for horizontal links. Note that, the flight efficiency costs for upward links are two times larger than downward links.

In Figure 11, links with high cost of noise impact are calculated. Referring to Figure 3, they are mainly distributed around densely populated areas.

D. Feasible path between OD pairs

To compare the relative importance of three types of link cost, in equation ( 9), the different combination of weighting parameters are examined. The weighting parameters w h , w f , w n are set to be [1/3, 1/3, 1/3], [2/3, 1/6, 1/6], [1/6, 2/3, 1/6], and [1/6, 1/6, 2/3]. In these four cases, the first parameter setting treats all link costs as equivalent. Other cases consider that one of the link costs is more important than the other link costs.

As the previous steps can be performed offline, we focus on the computation time of feasible path generation. These four route networks cost 11.19s, 11.74s, 13.39s and 11.31s, respectively. In addition, the average cost for k-shortest paths between OD pairs is shown in Figure 12. The coordinates indicate the index of vertiports and delivery sites. For OD pairs involving some vertiports, such as the vertiports 0, 1, 2, 3, the average path cost for different weighting parameters remains low. The reason is that they are more accessible to most delivery sites through the UAM route network. On the contrary, the average path costs for OD pairs that involve vertiport 4, 5, 6 are always at a high level. In Figure 12b, the link cost for crossing hazardous airspace is considered more important. In this case, the peak value of average path cost is the lowest among the four cases, because most hazardous airspace is covered by aerodromes, which is less penalized. Figure 12c corresponds to the case in which the link cost of flight efficiency is more important. The pattern of average path cost in this case is similar to the first case (Figure 12a), but the peak value is much higher. This fact indicates that feasible paths between some OD pairs are less efficient, such as OD pairs in the lower left part of Figure 12c. In the last case, it can be seen from Figure 12d that paths connecting vertiport 4 have significant costs. The reason lies in the fact that vertiport 4 is located in a densely populated downtown area.

To visualize the shortest path in the network representation, Figure 13 illustrate the number of times that each link appears in the k-shortest paths given different weighting parameters. A link is removed from the network if none of the feasible paths pass through it. At a quick glance, all altitude layers in these four networks are utilized. Vertical links, as the connection of adjacent altitude layers, appear most in the route network. According to different weighting parameters, the route network is created by minimizing the criteria formulized in equation [START_REF]Alos global digital surface model[END_REF]. Slight differences including link counts and network topology can be observed in these networks. To sum up, the resulting UAM route network provides many travel alternatives, which can be used by air traffic flow management and air traffic assignment to mitigate the complexity and congestion of airspace.

CONCLUSION AND PERSPECTIVES

The methodology proposed in this paper aims to define a set of alternative routes in a UAM context. Based on various data and grid networks, relevant constraints and costs are considered to search for k-shortest paths with diversity. The feasibility of this approach is demonstrated by its application to a parcel delivery service for Singapore.

In future work, it would be interesting to estimate the density of operations and analyze the impact of bottlenecks using different network structures. Another important extension of the work could consist in considering mixed traffic, for instance, automated air taxis sharing the airspace with parcel delivery UAVs. In addition, future studies could investigate the complete information processing chain by integrating the methodology proposed here with an air traffic assignment method using the sets of paths resulting from the solution to the KSPD problem. APPENDIX A. Exact expression of noise impact from a link to a polygon Assuming that individuals are uniformly distributed in a polygon and there is a uniform traffic flow of UAV on each link, then fS(a) = 1/|S| and fe(b) = 1/de, where |S| is the area of polygon and de is the length of line segment e. Le can be simplified as:

Le = 1 |S|de e S L dist(a,b) dadb (13)
As stated by [START_REF] Census | Geographic Areas Reference Manual[END_REF], many geographic regions including census tracts, block groups, and blocks are largely delineated by streets and are close to the rectangle in shape. To efficiently calculate the multiple integrals in Le, the subzone can be approximated by an AMABB [START_REF] Mu | Distance in spatial analysis: measurement, bias, and alternatives[END_REF][START_REF]Computation of the distance between a polygon and a point in spatial analysis[END_REF]. The route network is rotated by θ according to the reference point of AMABB, where θ is the angle by which the Minimum Bounding Box (MBB) rotates to become an AMABB.

After the rotation, the notation for calculating the noise metric of horizontal links and vertical links is given in Numerical integration method is used to solve these integrals.

Figure 1 :

 1 Figure 1: ALOS World 3D DSM data covering Singapore

  5km of aerodromes: The airspace within 5km of the aerodrome area. • Danger areas: The airspace within which activities dangerous to aircraft may exist at specified times. • Prohibited areas: The airspace within which aircraft is prohibited, usually due to security purposes. • Restricted areas: The airspace within which aircraft is restricted in accordance with certain specified conditions. • Protected areas: The reserved airspace within which aircraft is not allowed for operation. To sum up, prohibited and protected airspace are completely not allowed for UAM aircraft. All types of airspace are shown in Figure bounded by successive coordinates 2.

Figure 2 :

 2 Figure 2: Singapore existing area limits for UAM flights

Figure 3 :

 3 Figure 3: Total resident population of Singapore by subzone

Figure 4 :

 4 Figure 4: Geofence with keep-out distance of 30m

Figure 5 :

 5 Figure 5: Calculation of minimum delivery area and grid-based discretization

Figure 6 :

 6 Figure 6: The distribution of obstacles with geofence in delivery area at altitude 150ft, 175ft, 200ft, 225ft, and 250ft

Figure 7 :

 7 Figure 7: Measurement of SPL from UAV to a location

Algorithm 1

 1 Greedy algorithm for KSPD problem Input: Air transport network: G, OD pairs: W, similarity function: sim(•, •), number of shortest path: k, and similarity threshold: τ p Output: k-shortest path for OD pairs: Ψ 1: procedure KSPD 2:

3 :PΨ

 3 ← Set of acyclic paths that connects OD pair w and satisfies the constraint (1-2) 4: p * ← arg min p∈ P C p 5: P ← P/{p * } 6: Ψ w ← {p * } 7: while |Ψ w | < k and | P| > 0 do 8: p * ← arg min p∈ P C p 9: P ← P/{p * } 10: if ∀p ∈ Ψ w , sim(p, p * ) < τ p then 11: Ψ w ← Ψ w ∪ {p * } ← {Ψ w |w ∈ W} 16: return Ψ 17: end procedure B. Route network The resulting multi-layer route network is shown in Figure 8. The vertiports, delivery sites, and waypoints are marked in the figure. The links in different altitude layers are marked in different colors. The resulting route network contains 5 altitude layers.

Figure 8 :

 8 Figure 8: The resulting multi-layer route network in Singapore urban airspace

Figure 9 :

 9 Figure 9: Normalized link cost for crossing hazardous airspace in route network

Figure 10 :

 10 Figure 10: Normalized link cost of flight efficiency in route network

Figure A. 1 .L h + 10 log 10 h 2 h 2 + 2 h 2 + 2 dudxdy ( 15 ) 2 h 2 +

 1222221522 We give the exact expression of noise impact from a link to the nearby subzones:If the link is vertical, which corresponds to Figure A.1b, (x-x 0 ) 2 +(y-y 0 ) 2 dhdxdy (14) If the link is horizontal (corresponds to Figure A.1a) and the xcoordinates of both endpoints are different, Le = max(x0,x1) min(x0,x1) xmax xmin ymax ymin Lh + 10log 10 h (x-u) 2 +(y-(u-x 0 )(y 1 -y 0 ) x 1 -x 0 -y0) If the link is horizontal and the x-coordinates of both endpoints are the same, Le = max(y 0 ,y 1 ) min(y 0 ,y 1 ) (x-x 0 ) 2 +(y-v) 2 dvdxdy (16) (a) Horizontal link (b) Vertical link

Figure A. 1 :

 1 Figure A.1: Notation for different type of links

TABLE I :

 I Basic topology features of different altitude layers and the whole route network Network Number of nodes Number of links Average degree Average clustering coefficient Degree of assortativity Normalized link cost of noise impact in route network

	Network layer at 150ft	75	216	5.7600	0.3138	-0.0004
	Network layer at 175ft	102	336	6.5882	0.3875	0.3613
	Network layer at 200ft	111	444	8.0000	0.4709	0.2139
	Network layer at 225ft	115	524	9.1130	0.5081	0.2775
	Network layer at 250ft	115	546	9.4957	0.5154	0.3383
	The whole route network 529	2218	8.3856	0.3789	0.3789
	Figure 11:

ACKNOWLEDGEMENT

This work is partially supported by the research project CONCORDE of the Defense Innovation Agency (AID) of the French Ministry of Defense (2019650090004707501).