UDC 551.465

GEOPHYSICS

E. N. MIKHAILOVA, A. I. FELZENBAUM, N. B. SHAPIRO

ON THE COMPUTATION OF STEADY SEA AND OCEAN CURRENTS

(Presented by Academician L. I. Sedov, 28 VII 1965)

Consider the problem of determining a steady current caused by wind and climatological factors in a sea or ocean basin. Since the horizontal dimensions of the basin considerably exceed its depths, we have a kind of boundary layer at a solid boundary—the bottom of the basin. Hence there arises the possibility, first, of adopting the condition of hydrostatics and, second, in the equations of horizontal motion, of restricting ourselves to taking into account only the vertical exchange of momentum.

Introducing the integral stream function \(\psi\) and eliminating the pressure \(p\) by means of the hydrostatic condition * \((^1)\), we write the original system of equations and the boundary conditions in the form

\[ A\frac{\partial^2 u}{\partial z^2}+\Omega v = -g\frac{\partial \zeta}{\partial x} +\frac{g}{\rho_0}\int_0^z \frac{\partial \rho}{\partial x}\,dz +u\frac{\partial u}{\partial x} +v\frac{\partial u}{\partial y} +w\frac{\partial u}{\partial z}, \]

\[ A\frac{\partial^2 v}{\partial z^2}-\Omega u = -g\frac{\partial \zeta}{\partial y} +\frac{g}{\rho_0}\int_0^z \frac{\partial \rho}{\partial y}\,dz +u\frac{\partial v}{\partial x} +v\frac{\partial v}{\partial y} +w\frac{\partial v}{\partial z}; \tag{1} \]

\[ s_x=-\partial\psi/\partial y,\qquad s_y=\partial\psi/\partial x; \tag{2} \]

\[ w=\int_z^H\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)\,dz; \tag{3} \]

\[ u\,\partial\rho/\partial x+v\,\partial\rho/\partial y+w\,\partial\rho/\partial z = \varkappa_1\partial^2\rho/\partial z^2+\varkappa_2\Delta\rho; \tag{4} \]

for \(z=\zeta\)

\[ A\,\partial u/\partial z=-T_x/\rho_0,\qquad A\,\partial v/\partial z=-T_y/\rho_0; \tag{5} \]

\[ a_\zeta\partial\rho/\partial z+b_\zeta\rho=\Gamma_\zeta; \tag{6} \]

for \(z=H\)

\[ u=v=0; \tag{7} \]

\[ a_H\partial\rho/\partial z+b_H\rho=\Gamma_H; \tag{8} \]

on the contour of the basin \(l\)

\[ \psi=F(l); \tag{9} \]

\[ a_l\partial\rho/\partial n+b_l\rho=\Gamma_l. \tag{10} \]

In the equations and boundary conditions (1)—(10), \(u, v, w\) are the components of the current velocity along the Cartesian coordinate axes (the \(X\) axis is directed eastward, \(Y\) northward, \(Z\) vertically downward; the origin of coordinates is on the undisturbed surface of the ocean); \(\Delta=\partial^2/\partial x^2+\partial^2/\partial y^2\) is the Laplace operator; \(g\) is the acceleration of gravity; \(\Omega\) is the Coriolis parameter; \(\rho\) is the density of seawater (\(\rho_0=\mathrm{const}\) is the mean density); \(\zeta\) is the level; \(H\) is the depth of the basin; \(l\) is the contour of the basin and \(n\) is the direction of the normal to it; \(A\) is the coefficient of vertical exchange of momentum; \(\varkappa_1\) and \(\varkappa_2\) are the coefficients of vertical and horizontal turbulent diffu-

* For simplicity we shall regard the atmospheric pressure as constant, thus neglecting the static level \((^1)\).

boundary; \(T_x\) and \(T_y\) are the components of the tangential wind stress;

\[ s_x=\int_{\zeta}^{H} u\,dz \quad\text{and}\quad s_y=\int_{0}^{H} v\,dz \]

are the components of the total transport; \(F(l)\) is a known function on the contour of the basin \(l\) (if the basin is closed, then \(F(l)=0\)). The functions \(a_\zeta, b_\zeta, T_\zeta, a_H, b_H, T_H, a_l, b_l, T_l\) are assumed to be known. It should be noted that, generally speaking, these functions cannot be prescribed arbitrarily, and in each particular case the question of their choice must be resolved separately.

Let us write equations (1) in complex form

\[ \frac{\partial^2 V}{\partial z^2}-j^2 V=G+L, \tag{11} \]

where \(V=u+iv\); \(G=-\dfrac{g}{A}(\partial \zeta/\partial x+i\,\partial \zeta/\partial y)\); \(j^2=i\Omega/A\); \(L=L_1+iL_2=L_{11}+L_{12}+i(L_{21}+L_{22})\), with

\[ L_{11}=\frac{g}{\rho_0 A}\int_0^z \frac{\partial \rho}{\partial x}\,dz, \qquad L_{21}=\frac{g}{\rho_0 A}\int_0^z \frac{\partial \rho}{\partial y}\,dz, \]

\[ L_{12}=\frac{1}{A}\left(u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}\right), \qquad L_{22}=\frac{1}{A}\left(u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}\right). \]

If the right-hand side of equation (11) is regarded as known, then the solution of this equation under boundary condition (5), approximately referred to the undisturbed sea surface, and condition (7) has the following form (we introduce the notation \(T=T_x+iT_y\)):

\[ V=\frac{T}{jA}\frac{\operatorname{sh} j(H-z)}{\operatorname{ch} jH} +\frac{G}{j^2}\left[\frac{\operatorname{ch} jz}{\operatorname{ch} jH}-1\right] +\int_0^z L(s)\frac{\operatorname{sh} j(z-s)}{j}\,ds - \]

\[ -\int_0^H L(s)\frac{\operatorname{sh} j(H-s)\operatorname{ch} jz}{j\,\operatorname{ch} jH}\,ds. \tag{12} \]

Integrating expression (12) over the limits from \(z=\zeta\) to \(z=H\) and neglecting the quantity \(\zeta\), which is small in comparison with \(H\), we obtain for the total transport in complex form \(s=s_x+is_y\) the expression

\[ s=\frac{T}{j^2A}\frac{\operatorname{ch} jH-1}{\operatorname{ch} jH} +\frac{G}{j^2}\left[\frac{\operatorname{th} jH}{j}-H\right] + \]

\[ +\int_0^H\int_0^z L(s)\frac{\operatorname{sh} j(z-s)}{j}\,dz\,ds -\frac{\operatorname{sh} jH}{j^2\operatorname{ch} jH} \int_0^H L(s)\operatorname{sh} j(H-s)\,ds. \tag{13} \]

Determining the real and imaginary parts in expressions (12), (13), we obtain

\[ u=NT_x+MT_y+\Theta\frac{\partial \zeta}{\partial x}+\Lambda\frac{\partial \zeta}{\partial y}+ \tag{14} \]

\[ +\frac{1}{2a}\int_0^z (\alpha_1\alpha_5-\alpha_2\alpha_6-\alpha_3\alpha_7+\alpha_4\alpha_8)\,ds +\frac{1}{2a}\int_0^H (-\alpha_9\alpha_5+\alpha_{10}\alpha_6+\alpha_3\alpha_7-\alpha_4\alpha_8)\,ds, \]

\[ v=-MT_x+NT_y-\Lambda\frac{\partial \zeta}{\partial x}+\Theta\frac{\partial \zeta}{\partial y}+ \]

\[ +\frac{1}{2a}\int_0^z (\alpha_1\alpha_6+\alpha_2\alpha_5-\alpha_3\alpha_8-\alpha_4\alpha_7)\,ds +\frac{1}{2a}\int_0^H (-\alpha_9\alpha_6-\alpha_{10}\alpha_5+\alpha_3\alpha_8+\alpha_4\alpha_7)\,ds; \]

\[ \begin{aligned} s_x&=nT_x+mT_y+\vartheta\,\partial\zeta/\partial x+\lambda\,\partial\zeta/\partial y+s_x^*,\\ s_y&=-mT_x+nT_y-\lambda\,\partial\zeta/\partial x+\vartheta\,\partial\zeta/\partial y+s_y^*, \end{aligned} \tag{15} \]

where

\[ s_x^*=\frac{1}{2a}\int_0^H\int_0^z (\alpha_1\alpha_5-\alpha_2\alpha_6-\alpha_3\alpha_7+\alpha_4\alpha_8)\,ds\,dz\;-\; \]

\[ -\frac{1}{2a^2}\int_0^H (\alpha_{11}\alpha_5-\alpha_{12}\alpha_6-\operatorname{ch}aH\sin aH\,\alpha_7-\operatorname{sh}aH\cos aH\,\alpha_8)\,ds, \]

\[ s_y^*=\frac{1}{2a}\int_0^H\int_0^z (\alpha_1\alpha_6+\alpha_2\alpha_5-\alpha_3\alpha_8-\alpha_4\alpha_7)\,ds\,dz\;-\; \]

\[ -\frac{1}{2a^2}\int_0^H (\alpha_{11}\alpha_6+\alpha_{12}\alpha_5-\operatorname{ch}aH\sin aH\,\alpha_8+\operatorname{sh}aH\cos aH\,\alpha_7)\,ds, \]

where

\[ \begin{aligned} \alpha_1&=\operatorname{sh}az\cos az+\operatorname{ch}az\sin az, & \alpha_2&=\operatorname{ch}az\sin az-\operatorname{sh}az\cos az,\\ \alpha_3&=\operatorname{ch}az\cos az+\operatorname{sh}az\sin az, & \alpha_4&=\operatorname{sh}az\sin az-\operatorname{ch}az\cos az,\\ \alpha_5&=L_1\operatorname{ch}as\cos as-L_2\operatorname{sh}as\sin as, & \alpha_6&=L_1\operatorname{sh}as\sin as+L_2\operatorname{ch}as\cos as,\\ \alpha_7&=L_1\operatorname{sh}as\cos as-L_2\operatorname{ch}as\sin as, & \alpha_8&=L_1\operatorname{ch}as\sin as+L_2\operatorname{sh}as\cos as,\\ \alpha_9&=r(\alpha_3\operatorname{sh}2aH-\alpha_4\sin 2aH), & \alpha_{10}&=r(\alpha_3\sin 2aH+\alpha_4\operatorname{sh}2aH),\\ \alpha_{11}&=(2r+1)\operatorname{sh}aH\sin aH, & \alpha_{12}&=(2r-1)\operatorname{ch}aH\cos aH. \end{aligned} \]

Expressions for the remaining coefficients are given in paper (¹).

Let us now solve relations (15) with respect to the slopes of the sea surface:

\[ \frac{\partial \zeta}{\partial x} =-m'T_x+n'T_y+\vartheta'(s_x-s_x^*)-\lambda'(s_y-s_y^*), \]

\[ \frac{\partial \zeta}{\partial y} =-n'T_x-m'T_y+\lambda'(s_x-s_x^*)+\vartheta'(s_y-s_y^*). \tag{16} \]

Eliminating the level \(\zeta\) from relations (16) by cross differentiation and taking formulas (2) into account, we obtain

\[ \vartheta'\Delta\psi+ (\partial\vartheta'/\partial x+\partial\lambda'/\partial y)\,\partial\psi/\partial x+ (\partial\vartheta'/\partial y-\partial\lambda'/\partial x)\,\partial\psi/\partial y = \]

\[ =\operatorname{div}\,[\,n'T+\lambda's^*\,]+\operatorname{rot}_z[\,m'T+\vartheta's^*\,]. \tag{17} \]

It is easy to see that in the case of the linear theory for a homogeneous fluid
\((L_{11}=L_{12}=L_{21}=L_{22}=0)\) we arrive at the known relations of the theory (¹).

Passing to the numerical method for solving the problem, let us, for simplicity, consider two cases: a homogeneous fluid with account taken of the nonlinearity of the equations of motion, and an inhomogeneous fluid without this account.

In the case of a homogeneous fluid \((L_{11}=L_{21}=0)\) the solution is found with the aid of relations (1), (2), (3), (5), (7), (9). First we solve equation (17) with boundary condition (9), putting \(s_x^*=s_y^*=0\). Then, by formulas (16), (14), (3), we compute \(\partial\zeta/\partial x\), \(\partial\zeta/\partial y\), \(u\), \(v\), \(w\). Next we find \(L_{12}\), \(L_{22}\), and, finally, \(s_x^*\) and \(s_y^*\). Substituting the obtained values of \(s_x^*\) and \(s_y^*\) into equation (17) and solving it, we obtain the second approximation, and so on. The computations are carried out until two successive systems of the sought functions differ from one another by less than a specified small number.

In the second case, when the nonlinear inertial terms are discarded in the right-hand side of equations (1) \((L_{12}=L_{22}=0)\) and at the same time the density \(\rho\) is regarded as variable, we must add to the relations used in the preceding case the density diffusion equation (4) and the boundary conditions (6), (8), (10). As in the first case, at first we put, for example, \(s_x^*=s_y^*=0\) and solve equation (17) with boundary condition (9). Then we compute \(\partial\zeta/\partial x\), \(\partial\zeta/\partial y\), \(u\), \(v\), \(w\). The obtained values \(u\), \(v\), \(w\) are substituted into equation (4) and it is solved

under the indicated boundary conditions. Having found the density \(\rho\), we first compute \(L_{11}\) and \(L_{21}\), and then \(s_x^*\) and \(s_y^*\). Substituting the result into equation (17) and solving this equation with the new right-hand side, we find the next approximation, and so on. It should be noted that, as a first approximation, one need not use the solution of the problem for a homogeneous fluid; for example, one may set \(u=v=w=0\), solve equation (4), find \(s_x^*, s_y^*\), substitute their values into the right-hand side of equation (17), and so forth. The general case, when both the nonlinearity of the equations of motion and the inhomogeneity of sea water are taken into account simultaneously, is considered analogously.

The equation for the integral stream function (17) is solved by relaxation, splitting, and sweep methods \({}^{(2)}\). The density-diffusion equation (4) is solved by the same methods. Instead of this equation one considers the equation

\[ \partial \rho / \partial \eta + u \, \partial \rho / \partial x + v \, \partial \rho / \partial y + w \, \partial \rho / \partial z = \varkappa_1 \partial^2 \rho / \partial z^2 + \varkappa_2 \Delta \rho, \tag{18} \]

where \(\Delta \eta\) is a relaxation parameter. The interval \([\eta,\eta+\Delta \eta]\) is divided into three equal parts, and on each of them one of the following three equations is solved:

\[ \begin{aligned} \tfrac{1}{3}\partial \rho / \partial \eta + u \, \partial \rho / \partial x &= \varkappa_2 \partial^2 \rho / \partial x^2,\\ \tfrac{1}{3}\partial \rho / \partial \eta + v \, \partial \rho / \partial y &= \varkappa_2 \partial^2 \rho / \partial y^2,\\ \tfrac{1}{3}\partial \rho / \partial \eta + w \, \partial \rho / \partial z &= \varkappa_1 \partial^2 \rho / \partial z^2. \end{aligned} \tag{19} \]

Each of equations (19) is solved by the sweep method. The solution of equation (18), upon reaching a stationary regime, gives the solution of equation (4) \({}^{(3,4)}\).

In conclusion, we note that the nonstationary problem can also be solved, provided only that in the equations of horizontal motion the terms \(\partial u/\partial t,\ \partial v/\partial t\) are neglected. In this case the problem reduces to integrating the system of equations for the level \(\zeta\) and for the density \(\rho\). A similar, but less general, problem—because of certain additional assumptions—is considered in the work of A. S. Sarkisyan \({}^{(5)}\), to whom the idea of jointly considering the equations for these functions belongs.

Marine Hydrophysical Institute
Academy of Sciences of the Ukrainian SSR

Received
26 VII 1965

REFERENCES

\({}^{1}\) A. I. Felzenbaum, Theoretical Foundations and Methods for Calculating Steady Sea Currents, Moscow, 1960.
\({}^{2}\) A. I. Felzenbaum, N. B. Shapiro, DAN, 168, No. 3 (1966).
\({}^{3}\) G. I. Marchuk, DAN, 155, No. 5 (1964).
\({}^{4}\) G. I. Marchuk, Numerical Methods for Solving Problems of Weather Forecasting and Climate Theory, Part 1, Novosibirsk, 1965.
\({}^{5}\) A. S. Sarkisyan, Foundations of the Theory and Calculation of Oceanic Currents, Leningrad, 1966.