Karl Osen

Accurate Conversion of Earth-Fixed Earth-Centered Coordinates to Geodetic Coordinates

A closed form algorithm for the exact transformation of Earth-Centered Earth-Fixed (ECEF) coordinates to geodetic coordinates is presented that is computationally fast, safe and accurate. Boosting the computational robustness of Jijie Zhu's algorithm, it reduces the worst case transformation error by up to 500 million times.

Zhu's Algorithm

In 1993 Jijie Zhu published [START_REF] Zhu | Exact Conversion of Earth-Centered Earth-Fixed Coordinates to Geodetic Coordinates[END_REF] the following closed form algorithm for the exact transformation of Earth-Centered Earth-Fixed (ECEF) coordinates to geodetic coordinates :

= (1) (2) (3) (4)
(5)

(6) (7) (8) (9) (10) (11)
(12)

(13) (14) (15)
Special case when w = 0:

(16) (17)
In 1994 Zhu showed [START_REF] Zhu | Conversion of Earth-Centered Earth-Fixed Coordinates to Geodetic Coordinates[END_REF] that equation (9) can be rewritten to have only one cube root operation, thereby speeding up the transformation:

(18) Where:

(19)

Operating Range

Zhu states [START_REF] Zhu | Exact Conversion of Earth-Centered Earth-Fixed Coordinates to Geodetic Coordinates[END_REF] that the transformation of Earth-Centered Earth-Fixed (ECEF) coordinates to geodetic coordinates is feasible provided the point does not fall within 43 km of the Earth's center, and that this minimum distance warrants that the following condition is met: (20) Provided condition (20) is met we find in (7) that since . Consequently the radicand in (8) is always positive. We may furthermore show that which implies that we won't take the roots of negative numbers in (9) nor divide by zero in (18).

The feasibility analysis of Zhu's algorithm is essentially a two-dimensional problem, and we therefore introduce the two-dimensional coordinate notations and .

Replacing m and n in (20) we find:

(21) Where:

(22)

Inequality ( 21) implies that the point should be outside the ellipsoid with equatorial radius A and polar radius B. Hence, Zhu's 43 km rule has a safety margin of 158.7 meters. The ellipsoid equation is:

(24) Ellipsoid (24) has the same eccentricity as the WGS84 ellipsoid, but the polar radius is the larger radius. Whereas Zhu's 43 km rule can be used for most applications, the fundamental operating range of Zhu's algorithm is related to the radius of curvature in the meridian [START_REF] Clynch | Radius of the Earth -Radii used in Geodesy[END_REF]. This radius is the curvature in the north-south direction of the surface of the WGS84 ellipsoid at , and is given by: (25)

RM has its minimum at equator and its maximum at the poles:

(26) (27)
The center of the circle corresponding to RM has the following coordinates:

(28) (29)
We have and .

The line L starts in and ends in . The length of L is RM. L is orthogonal to the surface of the WGS84 ellipsoid in and its slope relative to the equatorial plane is . The points and are located on opposite sides of the equatorial plane (unless ). Consider a superellipsoid with the same equatorial and polar radii as ellipsoid (24):

If

(30)

Any point satisfy equation (30), which implies that the starting point of any L is located on the surface of this superellipsoid.

From the derivatives of ( 28) and (29) we find:

(31) Equation (31) shows that L is tangent to the surface of the superellipsoid in . The line L consequently passes through the superellipsoid (unless ).

Zhu's algorithm transforms ECEF coordinates to geodetic coordinates so that and are located on the same side of the equatorial plane. Zhu's algorithm does not work if is located on two lines L (i.e. one in the northern hemisphere and one in the southern hemisphere), which is the case for all points inside the superellipsoid.

The points and are on the surface of the superellipsoid, but generate divide-by-zero operations in equations ( 12) and (18). A point should therefore be outside the superellipsoid for successful conversion:

(32) Equation (8) has a component that tells us directly if the transformation of a point is feasible:

(33)
Inequality ( 32) is linked to equation (33) as follows:

(34)

Replacing and q in (33) we find:

Replacing m and n in (35) gives:

(36)

The derivatives of H with respect to w and z are:

(37) (38)
Combining ( 30) and (36) shows that for all points on the superellipsoid. Equations ( 37) and (38) show that and for . For we have except for and when . For we have except for and when . This implies that equation (34) holds.

The following test can therefore be used to check if an ECEF point mathematically can be converted to geodetic coordinates using Zhu's algorithm:

(39)

Unfortunately condition (39) is insufficient to ensure proper computational behavior of Zhu's algorithm close to the singularities at and . To stay sufficiently far away from these singularities we can advantageously introduce an exclusion volume fully containing the superellipsoid (30). Tests show that any point outside a volume V defined by , having equatorial radius meters and polar radius meters, avoids all numerical problems and accuracy degradations related to the proximity of the superellipsoid (30).

(40)

Consequently the following test can be used in to check if an ECEF point can accurately be converted to geodetic coordinates using Zhu's algorithm:

(41)

V looks like a deformed ellipsoid squeezed in at mid-latitudes, having its smallest distance from the center of the Earth of meters at

. The volume of V is of an ellipsoid with equatorial radius and polar radius . The maximum negative convertible altitude due to the exclusion volume V is meters at degrees.

Avoiding Negative Radicands

Numerical noise in implementations of Zhu's algorithm may result in (extremely small) negative radicands. Taking the root of a negative number will cause a computer program to fail or abort. To avoid such failures one may consider using absolute value functions before a root operations.

In equation (8) noise induced negative radicands may occur at very small latitudes, and the equation should therefore be modified as follows:

(42)

In equation (10) noise induced negative radicands may occur at latitudes close to ±45.3 degrees, and the equation should therefore be modified as follows:

(43)

Fixing the Accuracy Problem

The accuracy of Zhu's equation is generally quite good, but degrades considerably at latitudes close to ±45.3 degrees [START_REF] Zhu | Conversion of Earth-Centered Earth-Fixed Coordinates to Geodetic Coordinates[END_REF]. This accuracy degradation occurs when m is close to n. Combining equations ( 3) and (4) we find: The intersections between the cones and the surface of the WGS84 ellipsoid occurs at , which corresponds to degrees.

The intersections between the cones and Zhu's 43 km operating limit sphere occurs at , which corresponds to degrees.

By combining equation ( 44) with the superellipsoid (30) we find that the intersections between the cones and the superellipsoid (30) occurs at , which corresponds to degrees.

The loss of accuracy is due to the accumulation of errors when calculating , which is the real root of this 4 th degree equation:

(45)
If is computed with unlimited accuracy then is obviously zero, but unfortunately this is rarely the case. However, since both and its derivative can be calculated one may use Newton-Raphson's method to find a correction that can be added to :

(46) (47)
The value of is well positive and varies very slowly provided condition (41) is met, thereby ensuring risk free division in equation (47) and excellent error reduction (the overall conversion error may in some cases be reduced more than 100 million times).

To fix the accuracy problem of Zhu's algorithm we may now rewrite equations (11) and (12) as follows:

(48) (49)

WGS84 Constants

To obtain optimal accuracy for coordinate conversions using 32-bit, 64-bit or 80-bit floating point arithmetic, the following constants can advantageously be used:

Enhanced Algorithm

The following closed form algorithm for the exact transformation of Earth-Centered Earth-Fixed (ECEF) coordinates to geodetic coordinates is computationally faster, safer and more accurate than Zhu's original algorithm [START_REF] Zhu | Exact Conversion of Earth-Centered Earth-Fixed Coordinates to Geodetic Coordinates[END_REF]: 1. For optimal performance pre-computed constants should be used, replacing divisions by multiplications whenever possible (e.g. divide-by-6 should be replaced by multiply-by-0.166666666…). 2. Feasibility test (73) can be reached using only multiplications, additions and subtractions.

Feasibility test (73) uses

. This value is required for accurate 32-bit floating point operation.

can be smaller for 64-bit and 80-bit floating point arithmetic.

4. The particular formulation in (74) ensures that C > 0 if H > 0 .

5. The arctan operators in ( 86) and (90) take two parameters to avoid divide-by-zero situations.

Test Results

This chapter presents the numerical performance of Zhu's algorithm and the Enhanced algorithm, and compares accuracies for different floating point formats.

Conversion accuracy is measured by transforming Earth-Centered Earth-Fixed (ECEF) coordinates to geodetic coordinates and back to ECEF coordinates . The conversion error is given by .

Measurements were made on an Intel Core i5-4200U CPU, with all coordinate transformation calculations done on the x87 FPU. Trigonometric operations used FSINCOS and FPATAN instructions, whereas root operations were handled with FSQRT, FYL2X and F2XM1 instructions. The same transformation code featuring the sequence (66) to (90) was used for 32-bit, 64-bit and 80-bit floating point formats, but with different FPU configurations. The performance of Zhu's algorithm was measured by setting in (83) and (84). The inverse transformation function used is described by the sequence (61) to (65).

Measurements were made with 10 billion randomly generated positions for the different altitude ranges and floating point formats, with latitudes up to ±90 degrees and longitudes up to ±180 degrees. Conversion errors are due to accumulation of computational errors during transformations, but also of rounding errors when storing into floating point destination variables. Each number must be rounded to the nearest representation offered by the floating point format. The worst case rounding error is half of the floating point resolution for a given number. The floating point resolution of a number is the absolute value change caused by flipping the least significant bit of the number's floating point significand.

The expected error at is the combined worst case rounding errors for latitude, longitude and altitude, expressed as a geometric distance. Suppose the resolution at is and transforms to and transforms to , then the expected error is given by . The expected maximum error permits distinguishing the performance of the transformation algorithms from the limitations of floating point resolution. The expected maximum error is found at latitude zero, longitude ±180 degrees, and at the highest applicable altitude. The most commonly used floating point format for applications featuring ECEF and geodetic coordinates is the 64-bit format. Even when using 64-bit floating point numbers for coordinate representation it is nevertheless advantageous to use 80-bit floating point numbers during coordinate transformation to maximize accuracy. This puts the maximum error of the Enhanced algorithm just above the expected maximum error at .

  the point is located on the line L. If the point is located on the outwards extension of L. The distance from to is .

  represents two cones originating in the Earth's center. One cone extends northwards and the other southwards, and they are symmetric around the z axis. The angle between the cone surfaces and the equatorial plane is degrees. Points that are very far from the Earth have m = n for degrees.

  constants are useful for transforming geodetic coordinates to ECEF coordinates:

100km 2000km 1.07308e-09 6.54519e-09 3.216 175.110 3.37770e+05 2000km 35000km 2.07138e-09 2.51457e-08 11.375 -153.468 3.47446e+07 35000km 37000km 3.65036e-09 2.56071e-08 -16.869 -153.367 3.60552e+07 350000km 410000km 3.14502e-08 2.16964e-07 43.478 -80.559 4.02887e+08 146000000km 153000000km 1.31535e-05 9.94751e-05 -8.960 113.231 1.51844e+11 Floating point bits: 80 Tests per altitude range: 10000000000 alt.min: alt.max: avg.err: max.err: lat: lon: alt: -6378km -1km 4.49966e-13 3.18313e-12 -1.091 -17.575 -6.27794e+06 -1km 15km 5.20708e-13 3.19134e-12 -31.764 -156.767 1.23675e+04 15km 100km 5.79222e-13 3.38063e-12 -1.112 -125.856 5.10078e+04 100km 2000km 5.83478e-13 3.15958e-12 1.689 161.752 1.37304e+06 2000km 35000km 1.55298e-12 1.31248e-11 -9.203 57.473 3.41968e+07 35000km 37000km 2.70145e-12 1.36121e-11 19.423 31.678 3.66195e+07 350000km 410000km 2.30795e-11 1.23490e-10 -3.595 -45.037 3.88163e+08 146000000km 153000000km 1.00831e-08 5.38559e-08 11.873 160.439 1.51738e+11 #define WGS84_A +6.37813700000000000000e+0006 /* a */ #define WGS84_INVF +2.98257223563000000000e+0002 /* 1/f */ #define WGS84_F +3.35281066474748071998e-0003 /* f */ #define WGS84_INVA +1.56785594288739799723e-0007 /* 1/a */ #define WGS84_INVAA +2.45817225764733181057e-0014 /* 1/(a^2) */ #define WGS84_B +6.35675231424517949745e+0006 /* b */ #define WGS84_C +5.21854008423385332406e+0005 /* c */ #define WGS84_E +8.18191908426214947083e-0002 /* e */ #define WGS84_EE +6.69437999014131705734e-0003 /* e^2 */ #define WGS84_EED2 +3.34718999507065852867e-0003 /* (e^2)/2 */ #define WGS84_EEEE +4.48147234524044602618e-0005 /* e^4 */ #define WGS84_EEEED4 +1.12036808631011150655e-0005 /* (e^4)/4 */ #define WGS84_AADC +7.79540464078689228919e+0007 /* (a^2)/c */ #define WGS84_BBDCC +1.48379031586596594555e+0002 /* (b^2)/(c^2) */ #define WGS84_P1MEE +9.93305620009858682943e-0001 /* 1-(e^2) */ #define WGS84_P1MEEDAA +2.44171631847341700642e-0014 /* (1-(e^2))/(a^2) */ #define WGS84_P1MEEDB +1.56259921876129741211e-0007 /* (1-(e^2))/b */ #define WGS84_HMIN +2.25010182030430273673e-0014 /* (e^12)/4 */ #define WGS84_INVCBRT2 +7.93700525984099737380e-0001 /* 1/(2^(1/3)) */ #define WGS84_INV3 +3.33333333333333333333e-0001 /* 1/3 */ #define WGS84_INV6 +1.66666666666666666667e-0001 /* 1/6 */ #define WGS84_D2R +1.74532925199432957691e-0002 /* pi/180 */ #define WGS84_R2D +5.72957795130823208766e+0001 /* 180/pi */

Appendix A: Proof Related to Equations (30) and (36)

This appendix proves the following statement: Combining (30) and (36) shows that for all points on the superellipsoid.

We have:

We define:

We now have:

From (A3) we find:

Inserting (A6) and (A7) in (A4):

(A8)

Simplifying:

(A9) (A10)

Result:

(A11)

Appendix B: Zhu Algorithm Measurements

This appendix contains measurements of Zhu's algorithm.


P S Petrov
email: petrov@poi.dvo.ru

X Antoine
email: xavier.antoine@univ-lorraine.fr

Pseudodifferential adiabatic mode parabolic equations in curvilinear coordinates and their numerical solution

Keywords: parabolic equation method, multiple-scale method, curvilinear coordinates, computational ocean acoustics PACS: 43.30.+m, PACS: 43.20.Bi

Pseudodifferential mode parabolic equations in curvilinear coordinates are derived from the horizontal refraction equations. An energy-conserving splitstep Padé method for their numerical solution is developed. Several examples are considered including the propagation of sound in a perfect wedge, a whispering gallery formation near circular isobath, and the breaking of the whispering gallery waveguide in the vicinity of an inflection point of the isobath.

Introduction

The Mode Parabolic Equations (MPEs) approach is a promising tool for the simulation of sound propagation in shallow-water waveguides. MPEs in ocean acoustics were first derived independently by Collins in [START_REF] Collins | The adiabatic mode parabolic equation[END_REF] and Trofimov in [START_REF] Trofimov | Narrow-angle parabolic equations of adiabatic singlemode propagation in a horizontally inhomogeneous shallow sea[END_REF] (see also [START_REF] Petrov | Mode parabolic equations for the modeling of sound propagation in 3d-varying shallow water waveguides[END_REF]). In their pioneering works, mode coupling was neglected, and the derived equations were valid under the adiabatic assumption. Later on, the theory of narrow-angle MPEs was also generalized to the case of resonant mode coupling [START_REF] Trofimov | Narrow-angle parabolic equations of adiabatic singlemode propagation in a horizontally inhomogeneous shallow sea[END_REF][START_REF] Trofimov | A mode parabolic equation method in the case of the resonant mode interaction[END_REF][START_REF] Abawi | The coupled mode parabolic equation[END_REF]. In many important cases, adiabatic narrowangle MPEs even admit analytical solutions based on group-theoretical techniques [START_REF] Petrov | An explicit analytical solution for sound propagation in a three-dimensional penetrable wedge with small apex angle[END_REF][START_REF] Petrov | Analytical Lie-algebraic solution of a 3d sound propagation problem in the ocean[END_REF]. Recently, it was shown that in waveguides where the horizontal refraction is relatively strong the aperture of narrow-angle MPEs is insufficient to take all horizontal rays into account [START_REF] Trofimov | A mode parabolic equation method in the case of the resonant mode interaction[END_REF]. In such cases, Wide-Angle MPEs (WAMPEs) must be used to lead to an accurate simulation of the three-dimensional propagation of sound.

In many physical problems, some large portions of the acoustical energy emitted by a source in a shallow sea can be trapped in certain parts of the computational domain. For example, in some recent works, it was shown that, under certain conditions, whispering gallery waves localized in the vicinity of a curvilinear isobath can be excited in a shallow-water waveguide with bowl-like bottom [START_REF] Katsnelson | Whispering gallery waves localized near circular isobaths in shallow water[END_REF] (e.g. in a bay or a lagoon). It is also known that the sound field can be focused by an underwater canyon or trapped between circular fronts of internal waves [START_REF] Lynch | Acoustic ducting, reflection, refraction, and dispersion by curved nonlinear internal waves in shallow water[END_REF]. For such cases, it is natural to derive and solve mode parabolic equations in curvilinear coordinates.

In the present study, we combine the aforementioned motivations and derive PseudoDifferential [START_REF] Taylor | Pseudodifferential operators[END_REF][START_REF] Hörmander | The analysis of linear partial differential operators[END_REF] MPEs (PDMPEs) in orthogonal curvilinear coordinates. These equations have an arbitrarily wide aperture in the horizontal plane. We also propose a numerical method for solving PDMPEs on unbounded domains in the horizontal plane. This requires an artificial truncation of the computational domain by using Perfectly Matched Layers (PMLs). The numerical solution consists in a marching scheme, where a propagator associated with the PDMPE is evaluated at each step. The propagator is in turn a pseudodifferential operator that is replaced in practice by its Padé approximation. This technique was developed for similar equations complemented by periodicity condition in [START_REF] Antoine | Computing high-frequency scattered fields by beam propagation methods: A prospective study[END_REF] in the context of a scattering problem.

The plan of the paper is the following. In Section 2, we introduce a modal representation of an acoustical field lying in a 3D shallow-water waveguide and provide the horizontal refraction equation. The pseudodifferential mode parabolic equation in curvilinear coordinates is given in Section 3. Its numerical solution, based on a marching algorithm, is detailed in Section 4. Three numerical applications are developed in Section 5. Finally, we conclude in Section 6.

Modal representation of an acoustical field in a 3D shallow-water waveguide and the horizontal refraction equation

In this Section, we recall some key notions and equations related to adiabatic theory of sound propagation in 3D shallow-water waveguides. The physical conditions for the use of adiabatic approximation, as well as associated shortcomings and implications, are considered in many studies, e.g. in [START_REF] Jensen | Computational ocean acoustics[END_REF][START_REF] Pierce | Extension of the method of normal modes to sound propagation in an almost-stratified medium[END_REF]. Therefore, we do not discuss them here.

MPEs are usually derived [START_REF] Collins | The adiabatic mode parabolic equation[END_REF][START_REF] Trofimov | Narrow-angle parabolic equations of adiabatic singlemode propagation in a horizontally inhomogeneous shallow sea[END_REF] from the three-dimensional Helmholtz equation for a sound field in shallow-water produced by a monotone point source located at x = x s , y = y s , z = z s , i.e.

∂ 2 P ∂x 2 + ∂ 2 P ∂y 2 + ∂ 2 P ∂z 2 + ω 2 c 2 P = -δ(x -x s , y s , z -z s ) , (2.1)
where ω = 2πf is the cyclic frequency and c denotes the sound speed. In order to properly formulate the sound propagation problem, some interface and boundary conditions are required to be set up in addition to the Helmholtz equation (in particular, in shallow-water the continuity conditions for sound pressure and particle velocity at the water-bottom interface are of great importance). The complete formulation can be found in many textbooks and articles [START_REF] Petrov | Mode parabolic equations for the modeling of sound propagation in 3d-varying shallow water waveguides[END_REF][START_REF] Petrov | An explicit analytical solution for sound propagation in a three-dimensional penetrable wedge with small apex angle[END_REF][START_REF] Jensen | Computational ocean acoustics[END_REF].

The solution of (2.1) can be represented in the form of the modal expansion

P (x, y, z) = Nm j=1 A j (x, y)φ j (z, x, y) , (2.2)
where A j (x, y) are modal amplitudes, and φ j (z, x, y) are (local) mode functions computed for the media parameters (including the sound speed profile) at the point (x, y). Mode functions can be found from acoustical spectral problems [START_REF] Petrov | An explicit analytical solution for sound propagation in a three-dimensional penetrable wedge with small apex angle[END_REF][START_REF] Jensen | Computational ocean acoustics[END_REF][START_REF] Katsnelson | Fundamentals of shallow water acoustics[END_REF]. Numerous methods for its solution have been developed in the literature [START_REF] Jensen | Computational ocean acoustics[END_REF].

Under the assumption on the adiabaticity of the sound propagation, it can be shown that modal amplitudes satisfy the so-called Horizontal Refraction Equations (HREs) [START_REF] Collins | The adiabatic mode parabolic equation[END_REF][START_REF] Jensen | Computational ocean acoustics[END_REF]]

∂ 2 A j ∂x 2 + ∂ 2 A j ∂y 2 + k 2 j A j = -δ(x)δ(y) φ j (z s ) ρ w , (2.3)
where k j = k j (x, y) is the horizontal wavenumber corresponding to the eigenfunction φ j (z, x, y) and ρ w is the density in the water column.

Pseudodifferential mode parabolic equation in curvilinear coordinates

In many cases, it is natural to derive and solve MPEs in curvilinear coordinates. For example, in [START_REF] Katsnelson | Whispering gallery waves localized near circular isobaths in shallow water[END_REF], the authors considered whispering gallery waves formed in the vicinity of a curvilinear isobath due to the horizontal refraction of acoustical waves by the sloping bottom. Such waves can be accurately described in an orthogonal coordinate system where one of the coordinate lines coincides with the isobath (e.g., if the isobath is circular then the polar coordinate system should be used). As mentioned in the introduction, there also exist other propagation problems of shallow-water acoustics where it is natural to use curvilinear coordinates (e.g., the ducting of sound between the curved fronts of internal waves). For the rest of this study, we will mostly associate our coordinates with curved isobaths to have a concrete example where our results can be applied. However, it should be clear for the reader that in the general case one family of the curvilinear coordinate lines must coincide with level sets of the function k j = k j (x, y). If sound speed in the water is assumed to be constant, these level sets coincide with isobaths. Furthermore, our approach is also applicable to the solution of PDMPEs in Cartesian coordinate systems that are sufficient to handle the majority of sound propagation problems of shallow-water acoustics.

Assume that the source and the receiver are located on the same circular isobath where whispering gallery waves are excited. Our ultimate goal is to simulate both these waves creeping along an arc of the circle and direct-path propagation from the source to the receiver along the circle's chord. In order to tackle this problem, as well as more general cases where isobaths have more complicated forms, we need a WAMPE with an aperture of ±45 • that is solved along the isobath.

Let us assume that the isobath that supports whispering gallery waves in its vicinity is given by a parametric equation

x = X(s) , y = Y (s) , (3.1)
where s is the arclength. We start by introducing an orthogonal coordinate system (s, ξ) in its vicinity, where the ξ-axis is aligned along the depth gradient ∂h ∂x , ∂h ∂y and has the opposite direction (clearly, the ξ-axis has the same direction as the unit normal vector n = (Y (s), -X (s)) to the isobath (3.1). Any point (s, ξ) in the vicinity of the isobath has the Cartesian coordinates

x = X(s) + nY (s) , y = Y (s) -nX (s) . (3.2)
In the curvilinear coordinates (s, ξ), the HRE (2.3) writes as

1 γ ∂ ∂ξ γ ∂A j ∂ξ + 1 γ ∂ ∂s 1 γ ∂A j ∂s + k 2 j A j = 0 , (3.3)
where we omitted the source term on the right-hand side (this it typical for the derivation of PEs of various types). In (3.3), the variable coefficient γ = γ(s, ξ) = 1 + ξκ(s) is associated with the geometry of isobath, and κ(s) is the curvature of the latter. It can be computed as the length of the derivative of the unit tangent vector t = t(s) to the isobath

κ(s) = dt ds = d 2 X ds 2 2 + d 2 Y ds 2 2 . (3.4)
Derivation of PDMPEs corresponding to (3.3) can be greatly simplified by getting rid of the divergence form with respect to the arclength variable s. This can be accomplished by the substitution A j (s, ξ) = √ γv(s, ξ) 1 . The equation for v = v(s, ξ) has the following form

v ss + γ 2 v ξ ξ + γ 2 k 2 + 1 4 (γ ξ ) 2 + 1 2 γ ss γ - 3 4
(γ s ) 2 γ 2 v = 0 . (3.5)
It can be also written in the operator form

v ss = -L(s)v , (3.6)
where

L(s) = ∂ ∂ξ (γ 2 ) ∂ ∂ξ + γ 2 k 2 + 1 4 (γ ξ ) 2 + 1 2 γ ss γ - 3 4
(γ s ) 2 γ 2
is the transverse operator with respect to the propagation direction along the isobath.

From (3.6), a pseudodifferential PE (PDPE) can be written in the most obvious way as

v s = i L(s)v .
However, from general considerations, it is known that a better option is offered by an energy-conserving counterpart [START_REF] Antoine | Computing high-frequency scattered fields by beam propagation methods: A prospective study[END_REF][START_REF] Godin | Reciprocity and energy conservation within the parabolic approximation[END_REF][START_REF] Brekhovskikh | Wave propagation in a range dependent waveguide[END_REF][START_REF] Lu | Improving beam propagation method for TM polarization[END_REF][START_REF] Antoine | A construction of beam propagation methods for optical waveguides[END_REF] of the latter equation. It can be obtained by introducing a new unknown ϕ(s, ξ) = 4 L(s)v(s, ξ), and it can be easily shown that a PDPE for ϕ has the same form

ϕ s = i L(s)ϕ . (3.7)
Eq. (3.7) is called hereafter PseudoDifferential MPE (PDMPE). This nonlocal equation can be solved by using some Padé approximations for the square-root operator L(s).

The general recipe for the computation of an acoustical field can be outlined as follows:

1. Compute the modal wavenumbers k j (x, y) by solving acoustical spectral problems for all (x, y) (in practice, this can be greatly simplified by using mode perturbation theory). 2. Identify some suitable curvilinear coordinates (s, ξ) in such a way that the coordinate lines ξ = const coincide with the level sets of the function k j (x, y) (in many cases, it can be Cartesian coordinate system (x, y) chosen in a standard way with x-axis connecting the source and the receiver). 3. For each mode number j, solve Eq. (3.7) and compute v(s, ξ) by the formula v(s, ξ) := ( L(s)) -1/4 ϕ(s, ξ). 4. Evaluate A j (x, y) by transforming the function A j (s, ξ) = γv(s, ξ) back to the Cartesian coordinates (x, y). 5. Use the modal expansion (2.2) to compute the acoustical pressure P (x, y, z).

Numerical solution of the PDMPE: a marching algorithm

Following [START_REF] Antoine | Computing high-frequency scattered fields by beam propagation methods: A prospective study[END_REF], we introduce a uniform grid {s n } with step size h (i.e., such that s n+1 = s n + h). Equation (3.7) can be formally solved by a marching scheme [START_REF] Lu | Some techniques for computing wave propagation in optical waveguides[END_REF]

ϕ n+1 = exp(ik 0 h 1 + X)ϕ n := P ϕ n , (4.1)
where k 0 is the reference wavenumber, L = k 2 0 (1 + X), and the operator X is evaluated at the mid-point s n+1/2 = s n + h/2. The exponential operator P is called the one-way propagator corresponding to (3.7). Note that the respective propagators for narrow-angle MPE in many cases can be evaluated globally (for arbitrarily large steps) by various operator disentanglement analytical methods (even in the case of variable coefficients of special form, see [START_REF] Petrov | Analytical Lie-algebraic solution of a 3d sound propagation problem in the ocean[END_REF][START_REF] Petrov | Asymptotic solution for the problem of sound propagation in a shallow sea with the bathymetry described by a parametric quadratic function[END_REF]). By contrast, the more complicated case of (3.7) can be handled only numerically.

4.1. Approximation of the propagator P in Eq. (4. [START_REF] Collins | The adiabatic mode parabolic equation[END_REF] In this study we also use the same θ-approximant for the propagator P as described in [START_REF] Lu | Some techniques for computing wave propagation in optical waveguides[END_REF]

ϕ n+1 = P ϕ n ≈ e ik 0 h (d 0 + p k=1 d k 1 + b k X )ϕ n = e ik 0 h (d 0 ϕ n + p k=1 d k w n,k ) , (4.2)
where w n,k is obtained by solving the following auxiliary PDE

(1 + b k X)w n,k = ϕ n . (4.3)

Domain truncation and the PMLs

In the majority of practical problems of underwater acoustics, there exist areas, where the sound field has to be computed, that have no physical boundaries in horizontal directions. Thus, the problem of artificial domain truncations naturally arises. It can be solved by introducing a fictitious boundary at which either an artificial (or absorbing) boundary condition [START_REF] Antoine | A review of transparent and artificial boundary condition techniques for linear and nonlinear Schrödinger equations[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF] is set up or a Perfectly Matched Layer (PML) [START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Collino | Perfectly matched absorbing layers for the paraxial equations[END_REF] is attached to the computational domain. In this study, we choose the latter approach.

Assume that the solution of Eq. (3.7) is sought on the domain (ξ, s) ∈ Ω = [-d, d] × [0, s max ] (a curvilinear stripe of width 2d containing the given isobath (3.1)). We increase the width of the stripe by from each side and consider an initial-value problem for Eq. (3.7) on the extended domain

Ω = [-d -, d + ] × [0, s max ],
where the operator L is replaced by its PML version

LPML = 1 1 + iσ(ξ) ∂ ∂ξ γ 2 1 + iσ(ξ) ∂ ∂ξ + γ 2 k 2 + 1 4 (γ ξ ) 2 + 1 2 γ ss γ - 3

Discretization of the operator X

In order to solve the operator equations (4.3), we must discretize the operator X. It can be done by many techniques, but for the purposes of this study it is sufficient to use a standard second-order finite-difference method on a uniform grid {ξ }, where ξ +1 = ξ + τ . Thus, each operator equation (4.3) reduces to a tridiagonal system of linear algebraic equations.

Initial conditions for PDMPE

In computational problems of ocean acoustics, initial conditions for PEs (usually called PE starters) are chosen in such a way that the solution of the respective Initial-Boundary Value Problem (IBVP) approximates the solution of the Helmholtz equation (2.1) with a point source on the right-hand side. Many starters were developed for parabolic equations of various types (see e.g. [START_REF] Jensen | Computational ocean acoustics[END_REF]) in past decades. In our work, we use the so-called Greene's starter

A j (ξ, 0) = φ j (z s ) 2 √ π 1.4467 -0.8402k 2 j,0 ξ 2 e - k 2 j,0 ξ 2 1.5256 . (4.4)
Let us remark that this starter is probably non-optimal for PDMPE, as it was designed for a rational-linear PE. A development of the starter consistent with (3.7) will be the subject of a future work.

Examples

Sound propagation in a coastal wedge

In the first example, we consider the sound propagation in a wedge-shaped waveguide (see Fig. 1) that is a simplified model of a near-coastal area of a shallow sea. This is a standard test scenario used for validation of various 3D computational models in ocean acoustics. Although for many approaches the case of wedge with penetrable bottom is particularly challenging, this is not the case of MPEs. Indeed, difficulties associated with bottom properties are all confined to the solution of acoustical spectral problem (i.e., to the procedure of computation of k j and φ j (z, ξ, s)). On the other hand, perfectly reflecting bottom yields stronger manifestation of the horizontal refraction effect, and therefore is more suitable for validation of MPE-based models. Clearly, in this example, it is natural to use rectilinear Cartesian coordinates, where κ = 0 and γ = 1. However, to the best of our knowledge, the wedge problem has never been solved until now using wide-angle MPEs or PDMPEs. For our example, we take the bottom slope angle α = 1 • , the sound frequency f = 100 Hz and the water depth at the source location h 0 = 58 m. The sound speed in the water is 1500 m/s. The IBVP for Eq. (3.7) with initial conditions (4.4) is solved on the domain Ω = {(ξ, s)/|ξ| ≤ 3 km, 0 ≤ s ≤ 14 km}. At the boundaries ξ = ±3 km, PMLs with thickness = 300 m are added to the computational domain. The solution is obtained on a uniform grid with step size h = 2 m (the discretization step ∆ξ for the operator X is 0.5 m).

The contour plot of the computed solution is presented in logarithmic units in Fig. 2 (left panel). It is validated by a direct comparison with the solution provided by the method of source images [START_REF] Tang | On the method of source images for the wedge problem solution in ocean acoustics: Some corrections and appendices[END_REF] (see right panel of Fig. 2). It can be concluded from Fig. 2 that the two solutions are almost identical, although some small discrepancies are visible near the source (no- tice the sidelobes produced by PDMPE). For a more detailed comparison, we plotted the two solutions at ξ = 0 in Fig. 3, and the agreement appears to be nearly perfect.

Whispering gallery waves in the vicinity of a circular isobath (constant

curvature case) In the second example, we consider the sound propagation in a shallowwater with circular isobaths depicted in Fig. 4 (depth gradient must have the same direction as the normal vector to the isobath). As shown in [START_REF] Katsnelson | Whispering gallery waves localized near circular isobaths in shallow water[END_REF], under certain conditions, horizontal refraction in such a waveguide can result in the formation of Whispering Gallery (WG) waves [START_REF] Babič | Asymptotic methods in short-wavelength diffraction theory[END_REF] localized in the vicinity of an isobath. The properties of such waves and their modal structure was studied in [START_REF] Katsnelson | Whispering gallery waves localized near circular isobaths in shallow water[END_REF] using a WKBJ method. In this work, we present the results of the direct computation of the acoustical field produced by a point source in a waveguide of this type by the numerical solution of the PDMPE (3.7). Let us introduce the polar coordinate system (r, θ) in the horizontal plane in such a way that its center coincides with the curvature center of the isobaths. Let us assume that the bathymetry is described by a linear function h(r) on the interval [r 1 , r 2 ] (the same as in [START_REF] Katsnelson | Whispering gallery waves localized near circular isobaths in shallow water[END_REF]), i.e.

h(r) = h(r) =        h 1 if r ≤ r 2 , h 1 + h 2 -h 1 r 2 -r 1 (r -r 1 ) if r ∈ [r 1 , r 2 ] , h 2 if r ≥ r 2 .
(5.1)

The parameters are chosen as follows: r 1 = 5.8 km, r 2 = 6 km, h 1 = 26 m, h = h 2 = 24 m. In this example (as well as in the next one), the water column and the bottom are isovelocity layers, where the values of the sound speed are c w = 1500 m/s and c b = 2000 m/s, respectively. The densities in these layers are ρ w = 1 g/cm 3 , ρ b = 1.8 g/cm 3 , respectively. We assume that an omnidirectional point source of frequency f = 120 Hz is located at z = z s = 10 m, r = r s = 5800 m and θ = θ s = 0. In this case, it is natural to use curvilinear coordinates as

ξ = r -r s , s = r s θ , (5.2)
since WG waves are expected to be formed in the vicinity of the isobath r = r s (see Fig. 4).

The considered shallow-water waveguide supports three vertical modes for all values of r and θ. As shown in [START_REF] Katsnelson | Whispering gallery waves localized near circular isobaths in shallow water[END_REF], the WG-waveguide is formed in this case only for the highest-order vertical mode, i.e., for j = 3. In turn, the mode 3 component of the field P 3 (r, θ, z) = A j (r, θ)φ j (z, r, θ) exhibits a modal structure in the horizontal plane due to the formation of WG-modes. The energy transmitted by these modes is localized near the isobath r = r s , and does not undergo cylindrical spreading. Thus, in the far-field, the contribution of WG modes becomes dominant as the divergent part of the field fades away.

The method for the computation of the WG-component of P 3 (r, θ, z) is described in [START_REF] Katsnelson | Whispering gallery waves localized near circular isobaths in shallow water[END_REF]. It is based on the WKBJ approximation for the angular wavenumbers of WG modes and semi-analytical computation of radial mode functions. The solution magnitude in decibels 20 log 10 (4πP 3 ) for our example is shown in Fig. 5 (right panel) at the depth z = z s (the magnitude at 1 m from the source as the reference level). We also computed P 3 (r, θ, z = z s ) by solving PDMPE in the curvilinear coordinates (s, θ). The result is presented in Fig. 5 (left panel), and the comparison with the WG-component of P 3 calculated by the method from [START_REF] Katsnelson | Whispering gallery waves localized near circular isobaths in shallow water[END_REF] along the arc r = r s is shown in Fig. 6.

For a better illustration of this result, we also converted the PDMPE solution to the Cartesian coordinates (see Fig. 7). It can be seen that a certain part of the acoustical energy generated by the source is ducted in the vicinity of the isobath r = r s due to the horizontal refraction over the bottom slope. Although all horizontal rays corresponding to the mode 3 become bent due to this phenomenon, some of them do not get trapped in the WG-waveguide and leave the circle r = r 2 . After the respective energy is radiated into the outer space, the remaining part is formed exclusively by the WG-waves (as shown in Fig. 6).

In this example, we have κ = κ(s) = 1/r s . In addition, the comparison in

Dissipation of a WG wave: inflection point of an isobath

In the final example of our study, we consider a problem for which no analytical solution is available (apparently it cannot be also obtained by the MPE methods in a Cartesian coordinate system). In this case, radiation of the energy ducted in a WG-waveguide into the free-space will take place due to the presence of an inflection point of the isobath. This is an underwater acoustics counterpart of the WG propagation problem where the boundary of the domain has an inflection point first considered by Popov and his coauthors [START_REF] Popov | Problem of whispering gallery waves in a neighborhood of a simple zero of the effective curvature of the boundary[END_REF][START_REF] Popov | Numerical solution of the problem on whispering gallery waves in a neighborhood of a simple zero of the effective curvature of the boundary[END_REF] (both analytically and numerically). Some asymptotics of the field in shadow zones as well as the behaviour of the solution at infinity was also studied by Babich and Smyshlyaev [START_REF] Babich | Scattering problem for the Schrödinger equation in the case of a potential linear in time and coordinate. I. Asymptotics in the shadow zone[END_REF][START_REF] Babich | The scattering problem for the Schrödinger equation with a potential linear in time and in space. II. correctness, smoothness, behavior of the solution at infinity[END_REF] (see also the recent article [START_REF] Hewett | Contour integral solutions of the parabolic wave equation[END_REF] and references therein).

It is known that usual whispering gallery waves are formed near the reflective boundary of a convex domain. If at some point the curvature of the boundary changes its sign (i.e., the adjoint part of the domain is no longer convex), then the waveguide for WG waves ceases to exist, and the lion's share of their energy is radiated into the interior of the domain in the form of a searchlight beam. At the same time, minor part of the ducted energy penetrates into the shadow zone along the boundary after the inflection point. Clearly, one can expect similar behaviour of the WG waves trapped near a curvilinear isobath in a shallow sea.

Let us assume that a point source of frequency f = 120 Hz is placed at the isobath defined by the following formula in a Cartesian coordinate system (see also Fig. 8)

x = X(t) = x 0 -t , y = Ỹ (t) = a(t -x 0 ) 3 + ax 3 0 , (5.3)
where x 0 = 4000 m and a = 10 -8 m -2 are fixed constants (and t is a nonnatural parametrization of this curve). First, we introduce the arclength parameter s defined by

s = t 0 X2 (t) + Ỹ 2 (t)dt ,
rewrite the curve equations (5.3) as x = X(s), y = Y (s), and compute κ(s) using Eq. (3.4). The curvature radius 1/κ(s) gradually increases from 5.6 km at x = x 0 (i.e., from s = 0) to +∞ at x = 0. At x = 0, the isobath has an inflection point where the curvature κ(s) changes its sign (see the dashed lines in Fig. 8). Let us introduce the curvilinear coordinates (s, ξ) in the vicinity of this isobath, where s is the aforementioned arclength parameter of the latter, and ξ is the distance to the isobath from a given point in its vicinity. In the vicinity of the isobath, the depth varies only in the transverse direction, i.e.

h(ξ) =        h 1 if ξ ≤ ξ 2 , h 1 + h 2 -h 1 ξ 2 ξ if ξ ∈ [0, ξ 2 ] , h 2 if ξ ≥ ξ 2 .
(5.4)

Obviously, the dependence h = h(ξ) is exactly the same as in the previous example (although in this case it cannot be related to polar coordinates as in Eq. (5.1)), but the geometry of isobaths is more complicated due to the presence of an inflection point at x = 0. All the other environment parameters (sound speeds, densities, depths, etc) also have the same values as in the previous example.

In this case we also restrict our attention to the mode 3 component of the total field produced by the point source, since for the lower modes horizontal refraction does not lead to the formation of a WG waveguide. The contour plot of P 3 (x, y, z = z s ) in logarithmic units is presented in Fig. 8. The bottom is sloping between the isobaths shown in Fig. 8 by dashed lines and flat outside this area.

From Fig. 8, it is clear that the acoustical energy associated with mode 3 is concentrated in the vicinity of the isobath h = 26 m for x > 0 (thus, a WG is formed similarly to the previous example). For x < 0, i.e. after the inflection point, the WG waves form a searchlight beam and are radiated in the negative direction of the x-and y-axes. The diffracted wave further propagating into the shadow zone is also visible.

Conclusion

In this paper, we derived PDMPEs in orthogonal curvilinear coordinates. The derived equations can be used for the simulation of one-way propagation of sound in shallow-water waveguides with complicated bottom environment inhomogeneities that lead to complicated geometry of horizontal rays and to manifestation of various sound propagation effects. Some of these effects, e.g. whispering gallery formation and breaking, cannot be accurately simulated by using even very wide-angle MPEs in standard Cartesian or polar (cylindrical) coordinate systems since the propagation direction can change by 180 • at certain distance from the source. Due to the same reason, such problems also cannot be handled by three-dimensional parabolic equations.

At the same time, the examples in Sections 5.1 and 5.3 cannot be properly handled by narrow-angle MPEs similar to the equations derived in [START_REF] Babič | Asymptotic methods in short-wavelength diffraction theory[END_REF] and [START_REF] Popov | Numerical solution of the problem on whispering gallery waves in a neighborhood of a simple zero of the effective curvature of the boundary[END_REF]. Thus, WAMPEs or PDMPEs (similar to the one used in [START_REF] Antoine | Computing high-frequency scattered fields by beam propagation methods: A prospective study[END_REF]) are required in such cases.

It should be stressed that PDMPEs considered in our work are obtained under the adiabatic assumption and do not take mode coupling effects into account. However, for many practical problems, the horizontal refraction of sound is more important than the mode interaction, especially when the sound localization near an isobath is expected (or, more generally, near the level set of the horizontal wavenumber k j = k j (x, y)).

4 (γ s ) 2 γ 2 ,

 42 for a smooth function σ(ξ) which is increasing in ξ on [d, d + ], decreasing on [d -, d], and such that σ(ξ) = 0 for ξ ∈ [-d, d]. Thus, the operator LPML coincides with L on [-d, d], and outside this interval the derivatives ∂ ∂ξ in numerical scheme implements the solution of an initial-boundary value problem for Eq. (3.7) on the domain Ω with homogeneous Dirichlet boundary conditions of the form ϕ| ξ=±d± = 0 at the boundaries ξ = ±d ± of the PMLs. Inside the domain Ω, this solution accurately approximates the solution of Eq. (3.7) on the unbounded domain -∞ < ξ < ∞, s ≥ 0 provided that σ(ξ) is chosen in a suitable way. In the paper, we fix σ to σ(ξ) = σ 0 (ξ -d) 3 / 3 (in the examples presented below σ 0 = 5, = 300 m).

Figure 1 :

 1 Figure 1: Coastal wedge: a simplified model of a shallow sea with sloping bottom. The variable s is the distance from the source in the direction along the isobath, ξ is aligned along the depth gradient.

Figure 2 :

 2 Figure 2: Acoustical field in a wedge-shaped waveguide computed by the PDMPE (3.7) (left) and using the source images method [26] (right). The solutions have noticeable discrepancies near the source (especially in the area of PDMPE sidelobes) but perfectly coincide in the far-field.

Figure 3 :

 3 Figure 3: Comparison of the solutions obtained.

Figure 4 :

 4 Figure 4: An area of a shallow sea with circular isobaths. The variable s is the distance from the source along the circular isobath, ξ is aligned along the depth gradient. The depth is constant for ξ > ξ 2 and for ξ < ξ 1

Figure 5 :

 5 Figure 5: Acoustical field of the third mode A 3 (s, ξ) in a waveguide with circular isobaths computed by the PDMPE (3.7) (left) and using the WKB approximation for whispering gallery modes (right). The solutions perfectly coincide at long range, where the wave field is dominated by whispering gallery waves.

Figure 6 :

 6 Figure 6: Comparison of the two solutions shown in Fig. 5 along the arc r = r s (i.e., ξ = 0).

Figure 7 :

 7 Figure 7: Acoustical field in a an area of a shallow sea with circular isobaths where whispering gallery waves are formed due to horizontal refraction.

Fig. 6

 6 Fig. 6 indicates that the PDMPE accurately takes the curvilinear geometry of the domain into account.

Figure 8 :

 8 Figure 8: Radiation of the energy trapped in a whispering gallery waveguide into freespace. The isobaths h = 24m and h = 26m are shown by dashed lines.

Hereafter, we omit the subscript j indicating the mode number, and denote the partial derivatives by the subscripts s, ξ, etc.

Acknowledgements

The reported study was accomplished during P.S. Petrov's visit to the Institut Elie Cartan de Lorraine (Université de Lorraine, France) supported by the French government in the framework of the Metchnikov 2019 program. P.S. Petrov was also partly supported by the Russian Foundation for Basic Research (RFBR) under the contracts No. 18-35-20081 mol a ved and No. 18-05-00057 a. X. Antoine was supported by the French National Research Agency project NABUCO, grant ANR-17-CE40-0025.


Zineb Moumen
email: 1moumen.zineb@usmba.ac.ma

Najiba El

Amrani El

Manuela Tvaronavičienė
email: manuela.tvaronaviciene@vgtu.lt

Abderrahim Lahrach
email: abderrahim.lahrach@usmba.ac.ma

WATER SECURITY AND SUSTAINABLE DEVELOPMENT *

Keywords: water security, sustainable development, Green Morocco Plan, SDG6, Moroccan Dam's policy Additional disciplines: ecology and environment, environmental engineering

In the realms of Climate change impact, economic and population growth, and pollution, Water security, and sustainable development are paramount, for addressing the food-water-energy nexus and mitigating the impact of the cited drawbacks. The presented paper in the first part, emphasizes water security and sustainable development notions. In 2015 s the 2030 agenda were approved by 193 members of the United Nation, it consists of 17 goals (SDG) for environmental sustainability, social inclusion, economic development, peace, justice, good governance, and partnership. The 6th target of SDG is dedicated to clean water and sanitation since it is considered as one of the central focuses of Sustainable Development Goals. The second part is focused on Moroccan's water management plans. This latter elaborate a legal framework for water use, the law n° 36-15 that comprises the protection of human health through the regulation of exploitation, and provision of sanction. The green Morocco Plan (PMV) which aims first to develop modern agriculture aligned with the world global market requirement, and second to eliminate poverty by helping smallholder farming. The 140 dam's that Morocco counts provide strong support to economic growth. However, under the actual uncontrollable changing rate, Morocco is still threatened by low "water security", the access of drinkable water and sanitation services are improved but not totally covered, Moreover, the climate change

Introduction

From the starting of civilization evolving, water was not only a source of life but also a root of conflict and annihilation, because of extreme flood and drought. Through history, water trigger several wars, back to the world war II, the U.S bombard the North Vietnamese irrigation canals in the 1960s [START_REF] Jerome | A Brief History of Water Conflict[END_REF], in 1999s, the people protest against the privatization of water in Cochabamba, Bolivia [START_REF] Card | The Case of Privatization in Cochabamba, Bolivia Prepared for Julia Sagebien INTD 5003 Water Wars: The Case of Privatization in Cochabamba, Bolivia[END_REF], and the disobedience in Cape Town, because of insufficient water supply in 2012s [START_REF] Jerome | A Brief History of Water Conflict[END_REF]. In the other hand, the water excess may destroy humanity, a heavily unexpected rain may cause human and economic damages. For instance, in 1993s the Mississippi River overflew and inundated more than 9.3 million hectares, (Julie Bosman n.d.). According to Munich Re's natural catastrophe database (NatCatService), in Europe during the last 15 years, flash floods and storms have caused economic damages between 1 and 18 billion US dollars per incident. The Elbe and Danube floods in Germany and other parts of Eastern Europe in 2002 and 2013, caused more than 10 billion US dollars of damages ( Kottasova et al. n.d.), and according to a recent estimation, the risk of inundation in China is eight times higher than the global average [START_REF] Feng | Assessing the Inundation Risk Resulting from Extreme Water Levels under Sea-Level Rise: A Case Study of Rongcheng, China[END_REF].

According to the sixth report of intergovernmental panel of climate change (IPCC, 2014), a decrease of up to 20% in rainfall is predicted by the end of this century, and the increase in temperature is expected to reach 2.5 °C to 5.5 °C under the same scenario.

This paper introduces a broad definition of water security and sustainable development, the interaction between them, in the purpose of providing clear grasping of the concept. The second part is focused on Morocco as one of the countries in development that convoy sustainable development goals. It describes the legal framework of water, the role of the dam's policy and emphasizes the green Morocco plan (PMV). The last section points out the gaps in the cited strategies, and provide a recommendation that could be transposable to every country in development.

The hydrologic framework of water security

Only 2.5% of the total amount of water available of Earth is fresh water, 45 000Km3/year flow annually through the stream network to the sea [START_REF] Oki | Global Hydrological Cycles and World Water[END_REF].

Hydrological cycle assure the conservation of the water stock, the rainfall flows perennial and intermittent rivers, and supply groundwaters, the part that falls on the overland and vegetation forms the intercepted rain [START_REF] Perrin | Les Outils Du Projet : Modélisation Hydrologique et Représentation Des Usages[END_REF], the heat of sun energy, raises the water temperature, breaks the intermolecular links and transform it into a gas form [START_REF] Likens | Biogeochemistry of a Forested Ecosystem[END_REF].

The spatial distribution of precipitation is up to the geography and climate dynamics, there are regions with high water cyclicity, where extremes flood and droughts evoke risks of inundation and water scarcity.

303

The flow is the variation of water volume measured in cubic meter through the time, in another term, it is the circulation speed multiplied by the flowing area, this measure is relevant for the estimation of resident time of water molecule in a specific reservoir, by dividing its volume by the mean flux [START_REF] Oki | Global Hydrological Cycles and World Water[END_REF], in zones unaffected by anthropogenic activities, the residential time in water stream is around two and a half weeks (Institute of Industrial Science, Japan). In contrast, the recharge rate of groundwater is very slow according to the human time scale, the excessive uses lead to the exhausting of the resource [START_REF] Oki | Global Hydrological Cycles and World Water[END_REF].

Countries located in areas of high climate contrast (unpredictable floods and extensive drought) suffer from water insecurity. Grey and Sadoff (2007) emphasize the relation between hydrologic legacy and poverty rate, the study claims that countries which have an "easy" hydrologic legacy reach water security and have a high economic growth, in contrast, those with "difficult" hydrologic legacy are considered poor. This assumption is supported by [START_REF] Brown | Water and economic development: The role of variability and a framework for resilience[END_REF], who claim that high rainfall variability is statistically associated with lower capita incomes.

Certain aquifers or rivers could be shared between two or more countries, the world counts 263 transboundary rivers and 300 transboundary aquifers. For instance, Nile river is shared between Burundi, Central African, Republic, Egypt, Hala'ib triangle, Eritrea, Ethiopia, Kenya, Rwanda, Sudan, Abyei, South Sudan, United Republic of Tanzania, Uganda, Dem, and Republic of the Congo (Ntem and Melvin n.d.). The excessive use, the depletion, the pollution, the growing demand for water supply, imperil these shared resources and cause political and social conflict. 1948 Know 37 warfare about water (Thematic et al. n.d.), in the endeavor of resolving this encompassing issue, 295 international water agreement was signed (UNECE water convention, the legal framework for transboundary water cooperation ..) (Analysis n.d.).

Water security

Depending on the context and disciplinary perspectives on water use, the definition of "water security" that first emerged in the 1990s differs [START_REF] Cook | Water Security: Debating an Emerging Paradigm[END_REF], it is a coupled human-water system [START_REF] Srinivasan | A Dynamic Framework for Water Security[END_REF]) and a concept that provides a common language for all parties and components involved in water management (Arriens 2014), In the framework of UNESCO's International Hydrological Program's (IHP) Strategic Plan, Water Security was globally defined as the capacity of a population to safeguard sustainable access to adequate quantities of acceptable quality water for sustaining livelihoods, human well-being, and socio-economic development, for ensuring protection against water-borne pollution and water-related disasters, and for preserving ecosystems in a climate of peace and political stability [START_REF] Shrestha | Development of an Objective Water Security Index and Assessment of Its Association with Quality of Life in Urban Areas of Developing Countries[END_REF]. This implies that water security is a cycle with the implication of multiple interconnected and interdependent sectors or dimensions (hydrologic, geographic, economic, environmental, social, political, legal, financial…etc.) at local, national, regional, and global scales.

Water security concept is multidimensional, according to Van Beek and Arriens (2016), [START_REF] Falkenmark | The greatest water problem: the inability to link environmental security, water security and food security[END_REF] in his paper clarify the importance of bridging water security, with food and environmental security, as a major key to reach the sustainable development. Three key ones are to be consideredsocial equity, environmental sustainability, and economic efficiency (Table 1, Table 2).

INSIGHTS INTO REGIONAL DEVELOPMENT

Economic dimension Social dimension Environmental dimension

• increasing water productivity and conservation in all water-using sectors

• sharing economic, social, and environmental benefits in managing transboundary rivers, lakes, and aquifers.

• ensuring equitable access to water services and resources for all through robust policies and legal frameworks at all levels

• building resilience in communities in the face of extreme water events through both hard and soft measures.

• managing water sustainably as part of a green economy

• restoring ecosystem services in river basins to improve river health Source: [START_REF] Thematic | Water Security & the Global Water Agenda: A UN-Water Analytical Brief[END_REF] To assess water security at regional scale so many indexes were developed, some of them can be described in the table below; while these indicators are helpful, they are limited and don't take in consideration a lot of important dimensions. a national framework of water security measure applied to compare nations' water security performance in Asia and the Pacific region, The NWSI has five key dimensions: household, economic, urban, and environmental security and resilience to water-related disaster.

Asian development (ADB, 2016).

WII water insecurity index

It was applied on a regional scale in India and consisted of six key dimensions: resource, access, consumption, capacity, environment, and climatic stress Aggarwal, Punhani, and Kher (2014)

The concept of water shortage was quantifying first by Falkenmark in 1989, by developing the water crowding index, using a threshold value of water scarcity level, the critic limit occurs when the index is low to meet people's water needs (<1700 m3 cap-1 yr-1) [START_REF] Falkenmark | Meeting water requirements of an expanding world population [a discussion[END_REF]. The crowding index was the stepping stone of the evolving of water stress indicators, [START_REF] Raskin | Water Futures: Assessment of Long-range Patterns and Prospects[END_REF], [START_REF] Vörösmarty | Global water resources: Vulnerability from climate change and population growth[END_REF], [START_REF] Oki | Global assessment of current water resources using total runoff-integrating pathways[END_REF], and [START_REF] Alcamo | Global estimation of water withdrawals and availability under current and "business as usual" conditions[END_REF] develop other water stress indicator. The critical ratio compares the annual amount of available fresh water to withdrawal water [START_REF] Alcamo | Global change and global scenarios of water use and availability: An application of Water-GAP 1.0, Rep. A[END_REF][START_REF] Alcamo | World water in 2025-global modeling and scenario analysis for the world commission on water for the 21st centuryKassel World Water Series 2 Kassel[END_REF], [START_REF] Klepper | WARiBaS, water assessment on a river basin scale; a computer program for calculating water demand and water satisfaction on a catchment basin level for global-scale analysis of water stress Report 402001009[END_REF] develop a "water satisfaction ratio', [START_REF] Douglas | Characterizing water scarcity in Africa at different scales[END_REF] worked on the 'relative water stress index". All the above index converges to the same result, water stress is high when the used ration is superior to 0.4.

However, these indicators did not include the interaction between human and ecosystem [START_REF] Cook | Water Security: Debating an Emerging Paradigm[END_REF]. 305

INSIGHTS
Water resource faces other pressures, such as mismanagement, depletion, population, economic growth, and climate change impact, thus, the evolving of indicators that express better the relation between human-waterecosystem is crucial for a better grasp and then management [START_REF] Gain | Measuring Global Water Security towards Sustainable Development Goals Measuring Global Water Security towards Sustainable Development Goals[END_REF].

Sustainable development

The concept of sustainable development was first introduced in 1987 in the Report of the World Commission on Environment and Development called "our common future" as a "development that meets the needs of the present without compromising the ability of future generations to meet their own needs" (UN 1987). Sustainable Development (SD) is usually described as a three-dimension concept; environmental, economic and social (Giddings, Hopwood, & O´Brien, 2002). For a world called "sustainable", this concept must be achieved in all these dimensions.

A sustainable environment can be attained by reserving, improving and valuing the environment and natural resources in the long term, maintaining the major ecological balances, on the risks and on the environmental The UN entirely distinguishes the importance of managing water for sustainable development, and the need for the water and water-using sectors to collaborate and wisely use the world's limited water resources.

INSIGHTS
The SDG 6 namely sustainable water and sanitation for all, is allocated to water and comprises 11 global indicators, to track the progress.

Morocco between water insecurity and management strategies

Morocco is a northwest African country, cover an area of 710,850 km², the agricultural land represents 25 %, 80% of the total surface of the country receive less than 250 mm/year of rain (Africa et al. 2014), 80% evaporated and 0.6% is intercepted or infiltrated, the remained part is used for industry, irrigation and human consumption (Fig. The historical analysis of the hydrologic series data showed that morocco recorded during the last decencies, more than 20 dry periods (Africa et al. 2014(Africa et al. ), between 1982(Africa et al. -1983 the water deficit in Ziziguir watershed reach minus 80 mm, In 1995 the deficit was almost minus 100 mm in Loukous.

2.) (AQUASTAT). According to 2017 census

The frequency analysis of the chronologic data point out that the risk of drought is increasing from 5 years in 40 years to 4 years over 10 (199610 ( -200610 ( ) (Africa et al. 2014) ) These droughts have a direct impact on water supply and water ued for irrigation, for instance, in Tangier, there was a restriction of nearly 50% between 1993 and 1995, the deficit in the water allocated for irrigation rise from 16% to 90% between 199490% between -199590% between (Africa et al. 2014)). The law 36-15 includes cooperation users and the public authorities, aims to protect human health through the regulation of exploitation, distribution and sale of drinkable waters, provision for sanctions and the regulation of activities that may pollute water resources (Dahir n°1-16-113 du 6 kaada 1437 portant promulgation de la Loi n°36-15 relative à l'eau. 2016).

INSIGHTS INTO REGIONAL DEVELOPMENT

In addition to water laws, the policy of dams is one of the most effective and decent management strategies that Morocco adopt, it includes solutions for a set of critical drawbacks, and gives strong support to the economic growth. Now, Morocco counts 140 large dams, with a global capacity higher than milliards cubic meters, 13 are used for water transfer with a capacity of 200m3 /s, length 1100 km and a volume of 2.5 billion m3 / year [START_REF] El | La Politique et Les Bienfaits Des Barrages Au Maroc[END_REF], see The policy has a significant result concerning the above beneficiary sectors, the drinkable water and sanitation services in both rural and urban areas have increased, more than 50% of the total population have access to drinkable water and safely managed to basic sanitation service, (Fig. 6. and Fig. 7.). In the 2030s it is predicted to provide drinkable water and good sanitation services to the total population, which aligned with the 6 th sustainable development goal.

Pillar II:

The targets of this approach are the elimination of poverty by improving the performance of smallholder farming (around 560,000 farmers), specifically those living in the marginalized area, and this through three social projects (Aksebi 2012); -Reconversion projects that aim to switch from cereals to higher value products and less sensitive to water deficit, such as olive trees, almond trees, carob trees, cacti… -Diversification project and augmentation of territory product.

-Intensification project is concentrated on the supervising of farmers and providing the decent in order to enhance the yield and raise the production rate. The National Irrigation water saving program (PNEEI) aims to reduce water stress since it is the main limiting factor to the improvement of agricultural productivity. This program consists of conversion to localized irrigation over an area of nearly 550,000 hectares; an average equipment rate of nearly 55,000 ha/year (Ministry of Agriculture and Maritime Fisheries 2016).

In addition to drip irrigation, Morocco adopt other strategies for dealing with water stress issues, for example -The rainwater harvesting method based on the collection of water rain using different technics [START_REF] Hasnaoui | Plan de La Collecte Des Eaux Pluviales Plan de l ' Exposé[END_REF]), this later is legally managed by the decrees number 2-97-224 of 24/10/1997. 313

Soil composting help protecting the environment and especially water from the toxic substances impact, improving soil fertility and then agricultural productivity, reduce plant's pelt, the liquid extracted from the composting is used in drip irrigation (LABIDI 2016) and (Azin khalid, INRA d'agadir).

Conclusions

The water, food, energy, and ecosystems complex are a system in which there is plenty of interactions.

Water is the source of life and energy, which are crucial for food production, food is the basic element that maintains the ecosystem equilibrium. In contemporary society, energy is used for the water treatment process, groundwater extraction, and food production…

Regarding the environmental problems that the world face, the United Nation define 17 sustainable. Water is included in all forms of development, since it is the central focuses of Sustainable Development Goals, as it is at the heart of sustainable development and is crucial for socio-economic development, energy, and food production, healthy ecosystems and human survival (UN 2018).

Morocco is one of the countries the most threatened by climate change impact (IPCC) and water stress scarcity [START_REF] Gain | Measuring Global Water Security towards Sustainable Development Goals Measuring Global Water Security towards Sustainable Development Goals[END_REF]). Dam's policy and green morocco plan help the country to overlap issues related to water supply, increase the food security, decrease intern and extern migration, avoid the population centralization in big Cities, reduce social problems such as ignorance and poverty, provide a good life quality, However, it still facing problems, since the stress is still increasing, also the drip irrigation does not assure a strong production growth [START_REF] Tanouti | Squaring the Circle: Agricultural Intensification vs[END_REF], [START_REF] Alonso | Comparative Assessment of Irrigation Systems' Performance: Case Study in the Tri Ff a Agricultural District, NE Morocco[END_REF]).

In addition, the pesticide used to eliminate plant pelts help the growth of the agriculture productivity, but in other hand is not aligned with the ecosystem protection notion, the excessive uses of such product, halt the biodiversity and has a harmful sequel on human health, such as the deployment of cancers. The agricultural policy should be oriented toward a biologic production which respects the compromise between human benefice and ecosystems requirement since humankind is a part of the ecosystems and not the master of it.

Also, Morocco should integrate all dimension of sustainability and combine the economic progress with environment protection, by involving public and semi-public authorities (National office of drinkable water), privates companies for the elaboration of technical studies, International development institution (World Bank) and the most important invest in the scientific research which help providing advanced tools in terms of water management. For instance, build a national database including all the necessary information about soil, land use and land cover, hydrology…, these data could be introduced in mathematical tools for modeling, the obtaining results serve the decision makers to define the appropriate management strategies and help ecosystems to build a strong resilience capacity against climate change and natural hazard risks.

(

  [START_REF] Vörösmarty | Global threats to human water security and river biodiversity[END_REF] work on biophysical indicators to assess water security, (Lautze and Manthrithilake 2012) add the human dimension and assess water security in 46 Asian countries,[START_REF] Gain | Measuring Global Water Security towards Sustainable Development Goals Measuring Global Water Security towards Sustainable Development Goals[END_REF]) figure out the gap in the listed studies by providing the first global-scale assessment of water security, by the virtue of multicriteria analysis, the study conceptualize the security term in 4 majors part, each one is assessed by specific index (Fig.1.), and gather them into one normalized unique index.

Figure 1 .

 1 Figure 1. Water security global indicators Source: Gain et al 2016

  , the population is 35 M. The total yield of water is estimated by 22 Md m 3 (Africa et al. 2014), world water problems are downscaled on Moroccan scale, the resource is rare especially in the South and the South-eastern because of arid climate conditions, unequally distributed in space and time, facing a growing demand and climate change impact (Africa et al. 2014).

Figure 2 .Figure 3 .

 23 Figure 2. A subdivision of the Moroccan area, B: percentage of water loses from the total rainfall, C: Moroccan water uses. Source A: (Africa et al. 2014), B C: AQUAQTAT

Figure 4 .Figure 5 .

 45 Figure 4. Operational date of large dams in Morocco Source: Aquastat

Figure 6 .Figure 7 .

 67 Figure 6. The variation of sanitation rate in Morocco from 200 to 2015 Source: Household data

Table 1 .

 1 Economic, Social, Environmental dimension of water security

ISSN 2669-0195 (online)

http://jssidoi.org/jesi/ 2019 Volume 1 Number 4 (December) http://doi.org/10.9770/IRD.2019.1.4(2) 304

Table 2 .

 2 Economic, Social, Environmental dimension of water security

	Index	Meaning	Framework	developers
	NWSI	the national water
		security	index
		(NWSI)

indivisible and balance the three dimensions of sustainable development: the economic, social and environmental

  INTO REGIONAL DEVELOPMENT ISSN 2669-0195 (online) http://jssidoi.org/jesi/ 2019 Volume 1 Number 4 (December) http://doi.org/10.9770/IRD.2019.1.4(2) 306 impacts. A sustainable society can be maintained if it satisfies human needs and meet a social goal, by encouraging the participation of all social groups in health, housing, consumption, education, employment, culture…. etc., and finally a sustainable economy that aims for a developing growth and economic efficiency through sustainable production and consumption patterns (UN 1987). In 2000 the United Nations started the first steps by setting 8 original Millennium Development Goals (MDGs) to be achieved by 2015, to eradicate extreme poverty and hunger, achieve universal primary education, promote gender equality and empower women, reduce child mortality, improve maternal health, combat HIV/AIDS, malaria and other diseases, ensure environmental sustainability, and develop a global partnership for development. In September 2015 The Sustainable Development Goals (SDGs) were approved by the 193 UN Member States and are also known as Agenda 2030 (Transforming our World: the 2030 Agenda for Sustainable Development). The framework consists of 17 goals for environmental sustainability, social inclusion, economic development, peace, justice, good governance and partnership, the main issues for the world population in the 21st century. Each goal has several targets that better define its aims. While the MDGs aimed to reduce the existing issues mostly in under/developing countries by half, the new set of goals were set to eliminate them and has all countries as target (Sustainability Science 2018), and according to the agenda are "integrated and

(UN 2015)

.

Since water resources are embedded in all forms of development, Water is considered as one of the central focuses of Sustainable Development Goals, as it is at the heart of sustainable development and is critical for socio-economic development, energy and food production, healthy ecosystems and human survival (UN 2018),

Table 3 .

 3 Dams' benefices in drinkable water, irrigation, energy, drought and flood mitigation.

	Drinkable water	Irrigation			Energy		Drought	and	flood
							mitigation
	Dams satisfy 66% of the	1500,000	hectares	are	Current	hydroelectric	the management of 160 sites
	drinkable and industrial	irrigated from dams.		productivity is around 2.750	vulnerable to flooding by
	agglomerations water needs,				GWh / year,		dam construction.
	this rate is predicted to reach
	80% in 2020.				Hydroelectric	power	The lamination rate is high
					represents 10% of total	(100% in Abdelmoumen
					electricity production and	dam, 90.5 for Youssef ben
					will have to reach 14% by	each fin)
					the year 2020
			Source: (El and Directeur 2016)

* This research was supported by the project, which has received funding from the European Union's Horizon 2020 research and innovation programme European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme Marie Sklodowska-Curie Research and Innovation Staff Exchanges ES H2020-MSCA-RISE-2014 CLUSDEVMED (2015-2019) Grant Agreement Number 645730730 This research was supported by the project, which has received funding from the European Union's Horizon 2020 research and innovation programme European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme Marie Sklodowska-Curie Research and Innovation Staff Exchanges ES H2020-MSCA-RISE-2014 CLUSDEVMED (2015-2019) Grant Agreement Number 645730730 ://orcid.org/0000-0003-3413-6240 ://orcid.org/0000-0001-5603-1306

 ISSN 2669-0195 (online)