Determination of Dynamic Ocean Topography Using the Minimum Energy State

Geostrophic balance represents the minimum energy state in a linear system with conservation of potential vorticity. On the base of that, a minimum energy state (MES) method is proposed to determine dynamic topography (η) and in turn the absolute geostrophic velocity from hydrographic data through finding η to minimize a functional, { } 2 2 2 ( , ) ( ) 2 2 / min x y x y x y R G H C B f dxdy η η η η η η   = + + − →   ∫∫ . Numerical approach leads to a set of well-posed linear algebraic equations of η at grid points. Feasibility and advantage of the MES method are illustrated using the three-dimensional (T, S) data of the NOAA National Centers for Environmental Information (NCEI) World Ocean Atlas 2013 version 2 for the North Atlantic Ocean (100W-6W, 7N-72N) on 1× 1 grids.


Abstract
Geostrophic balance represents the minimum energy state in a linear system with conservation of potential vorticity. On the base of that, a minimum energy state (MES) method is proposed to determine dynamic topography (η) and in turn the absolute geostrophic velocity from hydrographic data through finding η to minimize a functional, Numerical approach leads to a set of well-posed linear algebraic equations of η at grid points. Feasibility and advantage of the MES method are illustrated using the three-dimensional (T, S) data of the NOAA National Centers for Environmental Information (NCEI) World

Introduction
Determination of ocean general circulation from temperature (T) and salinity (S) observations using the thermal wind relation is a classical problem. In the Cartesian coordinates with (x, y) in the horizontal and z in the vertical coordinates, the horizontal and vertical velocities are represented by (u, v) and w. As mentioned in [1], the quantities T, S are relatively easy to measure, and in contrast to velocity observations, the climatological signal in the (T, S) fields is less contaminated by energetic smaller-scale motions induced by eddies and waves. However, the hydrographic data only determine the baroclinic geostrophic currents (i.e., relative velocity between two depths). A reference velocity (i.e., velocity at reference depth z ref ) (u ref , v ref ) still needs to be determined by the inverse methods locally and non-locally.
The local inverse methods use the hydrographic data nearby the water column such as the β -spiral method [2] [3], the Bernoulli method [4], and the P-vector method [5]- [10]. The non-local inverse methods use the hydrographic data for the whole area such as the box method [11,12]. It was pointed out in [13] that the β -spiral method and the box method, no matter how different in appearance, are based on the same order of dynamical sophistication and differ from implicit assumptions about the scales of oceanic variability and different definitions of the smooth field to which the dynamical model pertains. ] at the ocean surface is calculated by and in turn the absolute geostrophic velocity at depth z is calculated by Here, ρ is the in-situ density; ρ 0 = 1025 kg m -3 is the characteristic density; g is the gravitational acceleration; f is the Coriolis parameter; and the Boussinesq approximation is used. On the base of geostrophic and hydrographic balances, and the conservation of the potential vorticity, the mechanical energy (kinetic and available potential energies) reaches minimum state as the ocean circulation 26 Determination of Dynamic Ocean Topography Using the Minimum Energy State takes the geostrophic velocity. With this minimum energy state, a new methodology has been developed to determine the dynamic topography (η) and in turn the absolute geostrophic velocity from hydrographic data. The main objective of this study is to show the theoretical background and technical part of this method. The rest of paper is outlined as follows. Section 2 depicts the minimum energy state (MES). Section 3 presents the MES inverse method including basic principles, numerical approach, and well-posed linear algebraic equations. Section 4 shows an example. Section 5 gives the conclusions.

Conservation of Potential Vorticity
Let the background fluid be motionless with the background density ( ρ ) being horizontally uniform and vertically increasing with depth linearly where N is the buoyancy frequency. The pressure field for the motionless background fluid is given by The density and pressure can be decomposed into background ( , p ρ ) and perturbation ( ρ , p) The linearized basic equations with the Boussinesq approximation are given by where 0 ρ is the characteristic density. Conservation of potential vorticity (q) can be derived from (6a)-(6e), where is a linearized potential vorticity. Here, (4a) and (4c) are used.

Extreme Energy State
The mechanical energy (sum of the kinetic and available potential energy) per unit mass is defined by The mechanical energy in a volume Ω is represented by Incorporation of the potential vorticity conservation constraint (8) leads to the functional where (12) and is the Lagrange multiplier, which is a function of space. If it were a constant, the integral would merely extremize energy subject to a given integral of potential vorticity, and rearrangement of potential vorticity would leave the integral unaltered. Extremization of the integral J leads to the three Euler-Lagrange equations, Substitution of (12) into (13), (14a, b) leads to (15) .
Differentiation of (16) with respect to z and use of (15) leads to

Basic Principle
Since E reaches the extreme state as the velocity takes the geostrophic velocity (u g , v g ) subject to a given potential vorticity. With the geostrophic and hydrographic balances, the large-scale the potential vorticity is calculated from the density field [5]. It can also be proved that E takes minimum state. However, it is beyond the scope of this study since the new inverse method is only on the base of the extremization of E (hereafter, the minimization is used instead of extremization) and the main objective here is to show the techniques. Substitution of (1)-(3) and (9) where H = H(x, y) is the depth of the water column, R is the horizontal area of the water volume, and the parameters (B, C) are defined by Use of the Euler-Lagrangian equation of the functional (20) leads to an elliptic partial differential equation for η, .

Numerical Approach
Let the three axes (x, y, z) be discretized into generalized (non-rectangular in horizontal) grid ( for a uniform rectangular horizontal grid. The subscripts in K i,j in (23a) indicates non-uniform water depth in the region. In a staggered grid, the density ( ) is evaluated at the original grid point (i, j, k). The absolute and baroclinic geostrophic velocities are evaluated at the center of the grid volume ( Figure 1). The dynamic topography (η) is at the original grid point (i, j) at the surface (k = 1). The functional G in (20) is discretized by , = 1, 2, ... , ; = 1, 2, ..., ; = 1, 2, ...,

Well-Posed Linear System
The minimization (20) can be conducted by (26) Substitution of (25) into (26) leads to a set of L linear algebraic equations with each equation having nine unknowns since the dynamic topography at the grid point η i,j is the corner point of four neighboring cells with (ε i-1, j-1 , ε i, j-1 , ε i-1, j , ε i, j ) evaluated at their centers and eight surrounding points The set of L linear algebraic equations can be written in the matrix form, where D is the L×L coefficient matrix, The Alternating-Direction Implicit (ADI) method described in Appendix B can also be used to solve the linear equation (30).

Treatment of Lateral Boundary
Three-types of grid cells are defined for the region R containing land ( Figure 3): sea cell (with no land mass), land cell (with no water mass), and boundary cell (with land and water masses). The dynamic topography ( ) is set to zero for the land cells and for the land parts of the boundary cells; and calculated for the sea and boundary cells. The absolute geostrophic velocity is calculated only at the sea cells. For the boundary cells, it is reasonable to set (33) because the absolute geostrophic velocity is not calculated. Substitution of (33) into (24) gives , for the water part in the boundary cells.

Treatment of Islands
Let the ocean domain Ω contain N islands. Since the normal velocity along the boundary of island is zero, it is reasonable to assume that the dynamic topography for a given island-n is a constant (no horizontal variability) η island-n . Let the Island-n have L boundary cells (Figure 4) with the value of η island-n . The dynamic topography at the sea cell grid point η island-n of the boundary cell is determined by the minimization (26), i.e.,

Example
The example (North Atlantic Ocean) described below is only to show the technical aspects of the proposed method. It is not our intention to investigate/review the dynamics of the North Atlantic Ocean circulation, which has already been done by many authors for a long time period.

Hydrographic Data
The   The three dimensional density is calculated from the WOA climatological annual mean temperature and salinity data (T i,j,k , S i,j,k ) using the international thermodynamic equation of seawater -2010, which is downloaded from the website: http://unesdoc.unesco.org/images/0018/001881/188170e.p df. The ocean depth data H i,j is downloaded from the NECI 5-Minute Gridded Global Relief Data Collection (ETOPO5) at the website: https://www.ngdc.noaa.gov/mgg/fliers/93mgg01.html. The numerical calculation is conducted using the WOA where ϕ is the latitude; r E = 6,371 km, is the earth radius, i = 1, 2, …, 95; j = 1, 2, …, 66 (i.e., I =95, J = 66). The parameters (A i,j , B i,j , C i,j ) from the data using (A3)-(A5). Then, the elements of the coefficient matrix D and the constant vector s are calculated from the values of (A i,j , B i,j , C i,j ) using (28) and (29), respectively. Use of (32) gives the dynamic topography vector η [or field η(x i , y j )].

Dynamic Topography and Absolute Geostrophic Velocity
The inverted dynamic topography η (Figure 7a) and associated surface geostrophic velocity vectors (Figure 7b The absolute geostrophic velocity vectors at 500m depth ( Figure 8a) are quite similar to the surface geostrophic velocity (Figure 7b) with reduced intensity. The maximum speed is around 10 cm s -1 (in the Gulf Stream) and the North Atlantic Subtropical Gyre shrinks in its north-south extension with increasing depth from the surface to 500m depth, which agrees with Olbers et al. (1985) and Chu (1995). The flow patterns at deeper ocean (1000m, 3500m, and 4500m) are shown in Figure 8b, 8c, 8d. The dominant feature at 1000m depth is the existence of southward flowing western boundary current along the east coast of North America and associated cyclonic gyre in the Sargasso Sea. This gyre is noticeable at deep depths (3500m, 4500m) with a maximum swirl speed around 2 cm s -1 . An evident southeastward flow (~ 2 cm s -1 ) is identified near Bahamas (south of 20 o N, 60 o W -50 o W) carrying water into the South Atlantic Ocean. Same as with the P vector method [5] (Chu 1995), a cyclonic-anticyclonic pair is also found in the deep water at the eastern (20 o W-40 o W) tropical (10 o N-30 o N) North Atlantic Ocean. The cyclonic eddy is in the south, and the anticyclonic eddy is in the north. The maximum swirl speed is around 2 cm s -1 . This feature is indeed strikingly similar to the map [14] of the absolute flow at 2000m depth (reproduced in [15]). Furthermore, our computation shows that this di-pole structure becomes more evident as the depth increases.

36
Determination of Dynamic Ocean Topography Using the Minimum Energy State The vertical cross-section of the zonal absolute geostrophic velocity component at 57 o W (Figure 9) shows the similar features identified using the box model [11] [16], the β spiral method [1], and the P-vector method [5]: (a) a strong upper ocean (above 1000 m depth) eastward flow between 30 o N-42 o N with a maximum speed around 12 cm s -1 , (b) a weak westward flow (~2 cm s -1 ) in deep layer below 3000m, (c) a near surface westward flow south of 30 o N (southern branch of the North Atlantic Gyre) with a maximum speed of 5 cm s -1 , and (d) banded structure in the deep layers (z < 1500m), where the current direction alternates on a horizontal scale around 2000km, which is wider than the results reported in [16].

Conclusions
On the theoretical base of minimum energy state of the geostrophic velocity with given potential vorticity, a straightforward MES method is proposed to determine the dynamic topography (η) and in turn the absolute geostrophic velocity from hydrographic data. This leads to a new mathematical problem how to find η to minimize the functional ( , ) x y G η η . In numerical approach, a set of corresponding linear algebraic equations is derived and solved to get η at grid points. Identification of η has several advantages in comparison to the existing methods that determine the velocity (u ref , v ref ) at a reference level. First, number of unknown variables reduces from 2 [i.e., (u ref , v ref )] to 1 [i.e., η]. Second, the geostrophic velocity is usually an order larger at the surface than at a deep reference level. Noise in (T, S) data has relatively larger impact on the deep reference-level velocity than on the surface geostrophic velocity. Third, the geostrophy is guaranteed for all the inverted three-dimensional absolute velocity using the MES method, but is not always guaranteed using the existing methods since it is difficult to prove that the reference velocity (u ref ,