RETRACTED ARTICLE: Ocean dynamic equations with the real gravity

Two different treatments in ocean dynamics are found between the gravity and pressure gradient force. Vertical component is 5–6 orders of magnitude larger than horizontal components for the pressure gradient force in large-scale motion, and for the gravity in any scale motion. The horizontal pressure gradient force is considered as a dominant force in oceanic motion from planetary to small scales. However, the horizontal gravity is omitted in oceanography completely. A non-dimensional C number (ratio between the horizontal gravity and the Coriolis force) is used to identify the importance of horizontal gravity in the ocean dynamics. Unexpectedly large C number with the global mean around 24 is obtained using the community datasets of the marine geoid height and ocean surface currents. New large-scale ocean dynamic equations with the real gravity are presented such as hydrostatic balance, geostrophic equilibrium, thermal wind, equipotential coordinate system, and vorticity equation.

Three forms of the gravity. Three forms of gravity exist in geodesy and oceanography: the real gravity g(λ, φ, z), the normal gravity [−g(φ)k], and the uniform gravity (−g 0 k, g 0 = 9.81 m/s 2 ). The real gravity (g) is a threedimensional vector field, which is decomposed into where g h is the horizontal gravity component; and (g z k) is the vertical gravity component. The normal gravity vector [−g(φ)k] (vertical vector) is associated with a mathematically modeled Earth (i.e., a rigid and geocentric ellipsoid) called the normal Earth. The normal Earth is a spheroid (i.e., an ellipsoid of revolution), has the same total mass and angular velocity as the Earth, and coincides its minor axis with the mean rotation of the Earth 7 . (1) (2) g = g h + g z k, g h = g i + g ϕ j www.nature.com/scientificreports/ Similar to g(λ, φ, z), the normal gravity [−g(φ)k] is the sum of the gravitational and centrifugal accelerations exerted on the water particle by the normal Earth. Its intensity g(φ) is called the normal gravity and determined analytically. For example, the World Geodetic System 1984 uses the Somiglina equation to represent g(φ) 8 where (a, b) are the equatorial and polar semi-axes; a is used for the Earth radius, R = a = 6.3781364 × 10 6 m; b = 6.3567523 × 10 6 m; e is the spheroid's eccentricity; g e = 9.780 m/s 2 , is the gravity at the equator; and g p = 9.832 m/ s 2 is the gravity at the poles. The uniform gravity (−g 0 k) is commonly used in oceanography.

OPEN
Real and normal gravity potential. Let (P, Q) be the Newtonian gravitational potential of the (real Earth, normal Earth) and P R (= � 2 r 2 cos 2 ϕ/2) ) be the potential of the Earth's rotation. Let V = P + P R be the gravity potential of the real Earth (associated with real gravity g) and E = Q + P R be the gravity potential of the normal Earth [associated with the normal gravity -g(φ)k]. The potential of the normal gravity [−g(φ)k] is given by The gravity disturbance is the difference between the real gravity g(λ, φ, z) and the normal gravity [−g(φ)k] at the same point 9 . The potential of the gravity disturbance (called the disturbing gravity potential) is given by Consequently, the centrifugal effect disappears and the disturbing gravity potential (T) can be considered a harmonic function. With the disturbing gravity potential T, the real gravity g (= g h + g z k) is represented by where ∇ h is the horizontal vector differential operator. The geoid height relative to the normal Earth (i.e., reference spheroid) is given by Bruns' formula 10 Equations (5), (6), (7) clearly show that the fluctuation of the marine geoid is independent of the Earth rotation and dependent on the disturbing gravity potential (T) evaluated at z = 0 only.
The disturbing static gravity potential (T) outside the Earth masses in the spherical coordinates with the spherical expansion is given by 11 where G = 6.674 × 10 −11 m 3 kg −1 s −2 , is the gravitational constant; M = 5.9736 × 10 24 kg, is the mass of the Earth; r is the radial distance with z = r−R; P l,m (sin ϕ) are the Legendre associated functions with (l, m) the degree and order of the harmonic expansion; (C l,m , C el l,m , S l,m ) are the harmonic geopotential coefficients (Stokes parameters with C el l,m belonging to the reference ellipsoid. From Eqs. (4) and (5) the potential of the real gravity is given by From Eq. (6) the real gravity is represented by The horizontal gravity component at the reference ellipsoid surface (z = 0) is obtained using (7) and (10) 12,13 An approximate 3D gravity field for oceanography. According to Eq. (8) (i.e., the spectral of the disturbing static gravity potential T), the ratio between T(λ, φ, z) to T(λ, φ, 0) through the water column can be roughly estimated by where H is the water depth. Since R is the radius of the Earth and more than 3 orders of magnitude larger than the water depth H. This leads to the first approximation that the surface disturbing gravity potential T(λ, φ, 0) is used for the whole water column, R r l C l,m − C el l,m cos m + S l,m sin m P l,m (sin ϕ), www.nature.com/scientificreports/ Since the deviation of the vertical component of the gravity (g z ) to a constant (−g 0 ) is 3-4 orders of magnitude smaller than g 0 , it leads to the second approximation With the two approximations, the near real gravity in the water column is given by Correspondingly, the potential of the real gravity is approximately given by where Eq. (9) is used.
Unexpectedly large horizontal gravity component. For simplicity, the local coordinates (x, y, z) are used from now on to replace the polar spherical coordinates with x representing longitude (eastward positive) and y representing latitude (northward positive). The local and spherical coordinate systems are connected by The EIGEN-6C4 model 14,15 , listed on the website http:// icgem. gfz-potsd am. de/ home, was developed jointly by the GFZ Potsdam and GRGS Toulouse up to degree and order 2190 to produce global static geoid height (N) dataset. Following the instruction, the author ran the EIGEN-6C4 model in 1° × 1° resolution for 17 s to get the global N ( Fig. 1) with mean value of 30.57 m, minimum value of -106.20 m, and maximum of 85.83 m. Two locations on the marine geoid are identified with N A = − 99.76 m at A (80° W, 3° N) in the Indian Ocean, and N B = 65.38 m at B (26° W, 45° N) in the North Atlantic Ocean, i.e., |ΔN| AB = 165.14 m. The big circle distance between A and B is about ΔL = 7920 km. Equation (11) shows that the corresponding |g h | at z = 0 is computed by With the geoid height data obtained from the EIGEN-6C4, N(x, y), the horizontal gravity components g x (Fig. 2a) and g y (Fig. 2b) at z = 0 are computed using Eq. (11). The magnitude of the horizontal gravity vector is calculated by .81 m/s 2 × 165.14 m 7.920 × 10 6 m = 20.45 mGal where f = 2Ωsinφ, is the Coriolis parameter; and U is the speed of the horizontal current. With Eq. (11) the C number at z = 0 is given by where ΔL is the horizontal scale of the motion.
The Ocean Surface Current Analysis Real-time (OSCAR) third degree resolution 5-day mean surface current vectors 17 on 26 February 2020 (Fig. 3a) was downloaded from the website https:// podaac-tools. jpl. nasa. gov/ drive/ files/ allDa ta/ oscar/ and the data at the same 1° × 1° grid points as the EIGEN-6C4 data were used. The data represent vertically averaged surface currents over the top 30 m of the upper ocean, which consist of a geostrophic component with a thermal wind adjustment using satellite sea surface height, and temperature and a wind-driven ageostrophic component using satellite surface winds. The histogram of the OSCAR current  (Fig. 4a). The histogram of (1/C 0 ) (Fig. 4b) indicates a positively skewed distribution with a long tail extending to values larger than 0.18. The statistical characteristics of (1/C 0 ) are 0.04238 as the mean, 0.04108 as the standard deviation, 1.447 www.nature.com/scientificreports/ as the skewness, and 4.481 as the kurtosis. Thus, the global mean of C 0 is about 23.6, which is coherent with the rough estimate given above using the geoid height data at the two points A and B (near 24).

New equations with the real gravity. Application of the Newton's second law of motion into the ocean
for the large-scale motions with the Boussinesq approximation leads to if the pressure gradient force, gravitation, and friction are the only real forces. Here, = � j cos ϕ + k sin ϕ , is the Earth rotation vector with Ω = 2π/(86,400 s) the Earth rotation rate; ρ is the density; ρ 0 = 1028 kg/m 3 , is the characteristic density; U = (u, v), is the horizontal velocity vector; ∇ is the three dimensional vector differential operator; and D/Dt is the total time rate of change; F = (F x , F y ), is the frictional force and usually represented using vertical eddy viscosity K,

The continuity equation is given by
where w is the vertical velocity. Substitution of (16) into (22) leads to where ω is the "vertical velocity" in the V-coordinate. Note that V is the static gravitational potential such that ∂V/∂t = 0.
Dynamic equations in the equipotential coordinate system. The horizontal pressure gradients can be computed from (27),

Substitution of (31) into (22) leads to the horizontal momentum equations
Extended geostrophic equilibrium and thermal wind relation. For steady state (no total derivative) without friction, Eqs. (32a) and (32b) become the extended equations for the geostrophic currents

Conclusion
The non-dimensional C number is used to identify relative importance of the horizontal gravity and the Coriolis force using the geoid height data (N) provided by the EIGEN-6C4 gravity model, and the Ocean Surface Current Analysis Real-time (OSCAR) surface current vectors on 26 February 2020. Unexpectedly large value of the C number (global mean around 24) may surprise both oceanographic and geodetic communities. It is very hard for my oceanographic colleagues to accept such strong horizontal gravity (at z = 0), which is an order of magnitude lager that the Coriolis force. Thus, in situ measurement of gravity may be needed for oceanography in addition to the routine hydrographic and current meter measurements. A geodesist may also be surprised that the horizontal gravity (z = 0) associated with the geoid height (N) varying between -106.20 m to 85.83 m (well accepted by the geodetic community), which was produced by a community gravity model (EIGEN-64C), generates the ocean currents with the intensity nearly 24 times as large as the currents identified from the ocean surface current analysis real-time (OSCAR) (oceanographic community product).
How to resolve such incoherency between the two communities on this issue becomes urgent. The new ocean dynamic equations including the horizontal gravity may provide a theoretical framework to resolve the incoherency since the potential of the real gravity (V) obtained from a geodetic gravity model is explicitly in the dynamical equations. Besides, use of the real gravity in the ocean dynamics may ultimately resolve some fundamental problems in oceanography such as reference level, and absolute geostrophic current calculation. A new dynamic system for large-scale oceanic motion with the real gravity is presented such as hydrostatic balance, geostrophic equilibrium, thermal wind, equipotential coordinate system, and the vorticity equation. Close collaboration between the oceanographic and geodetic communities helps the use of the real gravity in oceanography and the verification of the gravity model in geodesy with oceanographic data.

Appendix: Oblate spheroid coordinates versus polar spherical coordinates
This paper uses the polar spherical coordinates rather than the oblate spheroid coordinates for computational simplicity with a small error (0.17%) 18 . The oblate spheroid coordinates share the same longitude (λ) but different latitude (φ ob ) and radial coordinate (representing vertical) (r ob ). The relationship between the oblate spheroid coordinates (λ, φ ob , r ob ) and the polar spherical coordinates (λ, φ, r) is given by 18 where d is the half distance between the two foci of the ellipsoid. For the normal Earth, d = 521.854 km. The vector differential operator in the oblate spheroid coordinates is represented by www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. This is a U.S. Government work and not under copyright protection in the US; foreign copyright protection may apply 2021