21 Ocean circulation

Although this is a book on meteorology, there are three reasons to include a chapter on the ocean:

• Two thirds of Earth's surface are covered by the ocean. The ocean therefore constitutes a large part of the lower boundary condition of the atmosphere and is naturally crucial for the hydrological cycle.

• The ocean and atmosphere together form the climate system and interact with one another. Ocean circulation influences atmospheric climate by transporting heat and therefore modulates the position of fronts, cyclogenesis, and meridional climatological contrasts.

• The theoretical groundwork already developed for the atmosphere can be transferred rather easily to the ocean, so it would be a pity not to do so.

21.1 State variables of the ocean

The relevant state variables in the ocean are the three velocity components $u$, $v$, $w$, the surface elevation $\eta$, the thermodynamic variables $p$ and $T$, and the salt content (also called salinity) $S$. Alternatively, one may use quantities that can be derived bijectively from these, such as generalized velocities, the potential density, or the potential temperature. This is analogous to the state variables of the moist atmosphere without condensates, with salinity taking over the role played by humidity. One difference from the atmosphere is that the diabatic terms are very weak, especially in the deeper layers, owing to the weak radiation, the low flow speed and hence the low friction, and the absence of phase transitions. As a consequence, particles move essentially along isopycnals, i.e. surfaces of equal potential density.

21.2 Thermal equation of state

For liquids, in contrast to ideal gases, no simple, analytically derivable and sufficiently accurate thermal equations of state exist. In the case of seawater, the salt adds a further complication. A number of empirical equations of state for seawater have therefore been developed. The problem with empirical expressions, i.e. expressions fitted to measurement data, however, is their lack of thermodynamic consistency: suppose one has an empirical expression for the temperature as a function of volume, entropy and particle number, and another empirical expression for the pressure as a function of volume, entropy and particle number; then it is not guaranteed that the Maxwell relation Eq. (5.97) holds. To resolve this problem, the approach adopted for seawater is to use a formula for a thermodynamic potential and to compute other thermodynamic quantities by partial differentiation of this formula. This automatically ensures thermodynamic consistency. The modern standardized equation based on this principle is TEOS-10 (Thermodynamic Equation of Seawater 2010) [32].

This equation employs the Gibbs potential $G = G\left(T,p,m\right)$, or rather the specific Gibbs potential $g\left(T,p\right) = G\left(T,p,m\right)/m$. This transformation is possible because the Gibbs potential is an extensive state variable by virtue of Eq. (5.85). From this, the mass density can be computed by means of Eq. (5.95):

\[ \begin{align} \rho = \frac{m}{V} = \frac{1}{\frac{V}{m}} = \frac{1}{\frac{\partial g}{\partial p}} \end{align} \]

Here, the subscript $T,N$ or $T,m$, respectively, has been omitted, since these quantities are constant under this partial differentiation anyway, owing to $g = g\left(T,p\right)$.

21.3 Wind-driven circulation

The derivation of the Ekman spiral carried out in Sect. 17.3 can be transferred to a flat ocean floor without modification. At the ocean surface, however, the no-slip condition $\lim_{\mathbf{r} \to \partial V}\mathbf{v} = \mathbf{0}$ is absent; instead, it would be reasonable to assume that wind speed and current speed approach a common limiting value there. One would then have to solve the Ekman spiral simultaneously in the ocean and the atmosphere.

21.3.1 Sverdrup balance

There is, however, a simpler way to study the influence of wind on ocean circulation at a fundamental level. To this end, one first assumes the force balance given by Eqs. (17.35) - (17.36):

\[ \begin{align} f\mathbf{k}\times\mathbf{v}_{h} = -\nabla\phi + \frac{1}{\rho_0}\frac{\partial\mathbf{\tau}}{\partial z} \end{align} \]

It is now assumed that the vertical density gradients are so small that the factor $\frac{1}{\rho_0}$ can be absorbed into the partial derivative, i.e.

\[ \begin{align} f\mathbf{k}\times\mathbf{v}_{h} = -\nabla\phi + \frac{\partial}{\partial z}\left(\frac{\mathbf{\tau}}{\rho_0}\right). \end{align} \]

Applying the curl then yields

\[ \begin{align} \nabla\times f\mathbf{k}\times\mathbf{v}_{h} = \nabla\times\frac{\partial}{\partial z}\left(\frac{\mathbf{\tau}}{\rho_0}\right).\tag{21.4}\label{eq:wind-driven_circ_deriv_1} \end{align} \]

Using Eq. (B.53), this becomes

\[ \begin{align} \nabla\times f\mathbf{k}\times\mathbf{v}_{h} &= \left(\mathbf{v}_{h}\cdot\nabla\right)f\mathbf{k} - \mathbf{v}_{h}\left(\nabla\cdot f\mathbf{k}\right) + f\mathbf{k}\nabla\cdot\mathbf{v}_{h} - \left(f\mathbf{k}\cdot\nabla\right)\mathbf{v}_{h}. \end{align} \]

Projecting onto the vertical direction yields

\[ \begin{align} \mathbf{k}\cdot\nabla\times f\mathbf{k}\times\mathbf{v}_{h} &= \mathbf{k}\cdot\left(\mathbf{v}_{h}\cdot\nabla\right)f\mathbf{k} + f\nabla\cdot\mathbf{v}_{h}\nonumber\\ \Rightarrow \mathbf{k}\cdot\nabla\times f\mathbf{k}\times\mathbf{v}_{h} &= v\beta - f\frac{\partial w}{\partial z}. \end{align} \]

Integrating over the vertical interval $\left[z_B, z_T\right]$ yields

\[ \begin{align} \beta\int_{z_B}^{z_T}vdz - f\left(w_T - w_B\right) = \mathbf{k}\cdot\frac{1}{\rho_0}\nabla\times\left(\mathbf{\tau}_T - \mathbf{\tau}_B\right). \end{align} \]

With the definitions

\[ \begin{align} D \coloneqq z_T - z_B, & {} & \newoverline{v} \coloneqq \frac{1}{D}\int_{z_B}^{z_T}vdz \end{align} \]

this can be written more compactly as

\[ \begin{align} \beta D\newoverline{v} - f\left(w_T - w_B\right) = \mathbf{k}\cdot\frac{1}{\rho_0}\nabla\times\left(\mathbf{\tau}_T - \mathbf{\tau}_B\right). \end{align} \]

Since the entire water column is considered here, one can apply the simple boundary condition

\[ \begin{align} w_T = w_B = 0 \end{align} \]

This gives

\[ \begin{align} \beta D\newoverline{v} = \mathbf{k}\cdot\frac{1}{\rho_0}\nabla\times\left(\mathbf{\tau}_T - \mathbf{\tau}_B\right). \end{align} \]

If bottom friction is also neglected here, which is much smaller than wind stress because of the low flow speed in the deep ocean, one obtains the so-called Sverdrup balance, or Sverdrup relation:

In words, the Sverdrup balance states:“The vertically integrated meridional velocity equals the curl of the wind stress divided by the planetary vorticity gradient $\beta$.”

21.3.1.1 Stommel model

To parametrize the friction term at the bottom, one uses the Stommel model, which reads

\[ \begin{align} \mathbf{\tau}_B = r\newoverline{\mathbf{v}_{h}} \end{align} \]

with a constant $r > 0$. Since $w = 0$, one also has $\nabla\cdot \mathbf{v}_{h} = 0$, which is why the vertically averaged flow can be expressed in terms of a streamfunction $\psi$:

\[ \begin{align} \newoverline{u} = -\frac{\partial\psi}{\partial y}, & {} & \newoverline{v} = \frac{\partial\psi}{\partial x}\tag{21.14}\label{eq:stommel_stream} \end{align} \]

Thus one can write

\[ \begin{align} \beta D\frac{\partial\psi}{\partial x} + \mathbf{k}\cdot\frac{1}{\rho_0}\nabla\times\mathbf{\tau}_B &= \mathbf{k}\cdot\frac{1}{\rho_0}\nabla\times\mathbf{\tau}_T\nonumber\\ \Leftrightarrow \beta D\frac{\partial\psi}{\partial x} + \frac{r}{\rho_0}\Delta\psi &= \mathbf{k}\cdot\frac{1}{\rho_0}\nabla\times\mathbf{\tau}_T\nonumber \end{align} \]

with

\[ \begin{align} \epsilon \coloneqq \frac{\rho_0\beta D}{r}. \end{align} \]

Eq. (21.15) is the equation of motion of the Stommel model. From this, the ocean current can be diagnosed for a given wind forcing. To this end, Eq. (21.15) is solved on the set $\left[0, L_x\right] \times \left[0, L_y\right]$ with $L_x,L_y > 0$, i.e. for a rectangular basin. One starts from a wind stress of the form

\[ \begin{align} \tau_T^{(x)} &= -C\rho_aU^2\cos\left(\pi\frac{y}{L_y}\right),\\ \tau_T^{(y)} &= 0 \end{align} \]

where $\rho_a$ denotes the density of the atmosphere. From this it follows that

\[ \begin{align} \frac{1}{r}\mathbf{k}\cdot\nabla\times\mathbf{\tau}_T = -C\frac{\pi}{rL_y}\rho_aU^2\sin\left(\pi\frac{y}{L_y}\right). \end{align} \]

Inserting this into Eq. (21.15) yields

\[ \begin{align} \Delta\psi + \epsilon\frac{\partial\psi}{\partial x} &= -C\frac{\pi}{rL_y}\rho_aU^2\sin\left(\pi\frac{y}{L_y}\right).\tag{21.20}\label{eq:stommel_dynamics_mod} \end{align} \]

The boundary conditions are

\[ \begin{align} u\left(0,y\right) = u\left(L_x,y\right) &= 0,\nonumber\\ v\left(x,0\right) = v\left(x,L_y\right) &= 0, \end{align} \]

since no water is to leave the basin. From this it follows that

\[ \begin{align} \psi\left(0,y\right) = \psi\left(L_x,y\right) = \psi\left(x,0\right) = \psi\left(x,L_y\right) = 0. \end{align} \]

Furthermore, one now assumes that the depth is homogeneous and that $f$ depends linearly on latitude, so that $\beta$ is homogeneous ($\beta-$plane). Eq. (21.20) is linear, i.e. sums of solutions are again solutions. Moreover, this equation is inhomogeneous. One first solves the homogeneous part

\[ \begin{align} \Delta\psi + \epsilon\frac{\partial\psi}{\partial x} &= 0. \end{align} \]

One makes the ansatz

\[ \begin{align} \psi\left(x,y\right) = \left[A\exp\left(\lambda_1x\right) + B\exp\left(\lambda_2x\right)\right]\sin\left(\frac{m\pi y}{L_y}\right). \end{align} \]

Insertion yields

\[ \begin{align} &\left[A\lambda_1^2\exp\left(\lambda_1x\right) + B\lambda_2^2\exp\left(\lambda_2x\right)\right]\sin\left(\frac{m\pi y}{L_y}\right) - \frac{m^2\pi^2}{L_y^2}\left[A\exp\left(\lambda_1x\right) + B\exp\left(\lambda_2x\right)\right]\sin\left(\frac{m\pi y}{L_y}\right)\nonumber\\ &+ \epsilon\left[A\lambda_1\exp\left(\lambda_1x\right) + B\lambda_2\exp\left(\lambda_2x\right)\right]\sin\left(\frac{m\pi y}{L_y}\right) = 0\nonumber\\ &\Leftrightarrow\left[A\lambda_1^2\exp\left(\lambda_1x\right) + B\lambda_2^2\exp\left(\lambda_2x\right)\right] - \frac{m^2\pi^2}{L_y^2}\left[A\exp\left(\lambda_1x\right) + B\exp\left(\lambda_2x\right)\right] + \epsilon\left[A\lambda_1\exp\left(\lambda_1x\right) + B\lambda_2\exp\left(\lambda_2x\right)\right] = 0 \end{align} \]

Since $\exp\left(\lambda_1x\right)$ and $\exp\left(\lambda_2x\right)$ are linearly independent of each other, the sums of their respective prefactors must each vanish, i.e.

\[ \begin{align} \lambda_1^2 + \epsilon\lambda_1 - \frac{m^2\pi^2}{L_y^2} &= 0,\\ \lambda_2^2 + \epsilon\lambda_2 - \frac{m^2\pi^2}{L_y^2} &= 0. \end{align} \]

From this, the quadratic formula yields

\[ \begin{align} \lambda_1 &= -\frac{\epsilon}{2} + \sqrt{\frac{\epsilon^2}{4} + \frac{m^2\pi^2}{L_y^2}},\\ \lambda_2 &= -\frac{\epsilon}{2} - \sqrt{\frac{\epsilon^2}{4} + \frac{m^2\pi^2}{L_y^2}}. \end{align} \]

One now also needs a particular solution of Eq. (21.20) including the inhomogeneous part. For this, one chooses

\[ \begin{align} \psi\left(x,y\right) = C\frac{L_y}{r\pi}\rho_aU^2\sin\left(\pi\frac{y}{L_y}\right). \end{align} \]

The general solution therefore reads

\[ \begin{align} \psi\left(x,y\right) = \left[A\exp\left(\lambda_1x\right) + B\exp\left(\lambda_2x\right)\right]\sin\left(\frac{m\pi y}{L_y}\right) + \frac{CL_y}{r\pi}\rho_aU^2\sin\left(\pi\frac{y}{L_y}\right). \end{align} \]

$A$, $B$, and $m$ can be derived from the boundary conditions. The boundary conditions $\psi\left(x,0\right) = \psi\left(x,L_y\right) = 0$ are satisfied. From $\psi\left(0,y\right) = \psi\left(L_x,y\right) = 0$ for all $0\leq y\leq L_y$, one first obtains $m=1$, since otherwise the two sine terms cannot cancel each other. Hence

\[ \begin{align} \psi\left(x,y\right) = \left[A\exp\left(\lambda_1x\right) + B\exp\left(\lambda_2x\right) + \frac{L_y}{r\pi}C\rho_aU^2\right]\sin\left(\pi\frac{y}{L_y}\right).\tag{21.32}\label{eq:stommel_deriv_2} \end{align} \]

With $\psi\left(0,y\right) = \psi\left(L_x,y\right) = 0$, one now obtains the following linear system of equations for $A$ and $B$:

\[ \begin{align} A + B &= -\frac{L_y}{r\pi}C\rho_aU^2,\tag{21.33}\label{eq:stommel_deriv_1}\\ A\exp\left(\lambda_1L_x\right) + B\exp\left(\lambda_2L_x\right) &= -\frac{L_y}{r\pi}C\rho_aU^2 \end{align} \]

Multiplying the first equation by $\exp\left(\lambda_1L_x\right)$ and subtracting the result from the second, one obtains

\[ \begin{align} B\left[\exp\left(\lambda_2L_x\right) - \exp\left(\lambda_1L_x\right)\right] &= \frac{L_y}{r\pi}C\rho_aU^2\left(\exp\left(\lambda_1L_x\right) - 1\right)\nonumber\\ \Leftrightarrow B&= \frac{L_y}{r\pi}C\rho_aU^2\frac{\exp\left(\lambda_1L_x\right) - 1}{\exp\left(\lambda_2L_x\right) - \exp\left(\lambda_1L_x\right)} \end{align} \]

Inserting this into Eq. (21.33) yields

\[ \begin{align} A &= -\frac{L_y}{r\pi}C\rho_aU^2 - B = -\frac{L_y}{r\pi}C\rho_aU^2 - \frac{L_y}{r\pi}C\rho_aU^2\frac{\exp\left(\lambda_1L_x\right) - 1}{\exp\left(\lambda_2L_x\right) - \exp\left(\lambda_1L_x\right)}\nonumber\\ &= -\frac{L_y}{r\pi}C\rho_aU^2\frac{\exp\left(\lambda_2L_x\right) - \exp\left(\lambda_1L_x\right)}{\exp\left(\lambda_2L_x\right) - \exp\left(\lambda_1L_x\right)} - \frac{L_y}{r\pi}C\rho_aU^2\frac{\exp\left(\lambda_1L_x\right) - 1}{\exp\left(\lambda_2L_x\right) - \exp\left(\lambda_1L_x\right)}\nonumber\\ &= -\frac{L_y}{r\pi}C\rho_aU^2\frac{\exp\left(\lambda_2L_x\right) - 1}{\exp\left(\lambda_2L_x\right) - \exp\left(\lambda_1L_x\right)}. \end{align} \]

Inserting the results for $A$ and $B$ into Eq. (21.32) and further defining

\[ \begin{align} \gamma \coloneqq \exp\left(\lambda_2L_x\right) - \exp\left(\lambda_1L_x\right), \end{align} \]

one obtains

\[ \begin{align} \psi\left(x,y\right) &= \frac{L_y}{r\pi}C\rho_aU^2\left[\frac{1 - \exp\left(\lambda_2L_x\right)}{\gamma}\exp\left(\lambda_1x\right) + \frac{\exp\left(\lambda_1L_x\right) - 1}{\gamma}\exp\left(\lambda_2x\right) + 1\right]\sin\left(\pi\frac{y}{L_y}\right).\tag{21.38}\label{eq:stommel_solution} \end{align} \]

Inserting this into Eq. (21.14) yields

\[ \begin{align} u\left(x,y\right) &= -\frac{C\rho_aU^2}{r}\left[\frac{1 - \exp\left(\lambda_2L_x\right)}{\gamma}\exp\left(\lambda_1x\right) + \frac{\exp\left(\lambda_1L_x\right) - 1}{\gamma}\exp\left(\lambda_2x\right) + 1\right]\cos\left(\pi\frac{y}{L_y}\right),\tag{21.39}\label{eq:stommel_result_u}\\ v\left(x,y\right) &= \frac{L_y}{r\pi}C\rho_aU^2\left[\lambda_1\frac{1 - \exp\left(\lambda_2L_x\right)}{\gamma}\exp\left(\lambda_1x\right) + \lambda_2\frac{\exp\left(\lambda_1L_x\right) - 1}{\gamma}\exp\left(\lambda_2x\right)\right]\sin\left(\pi\frac{y}{L_y}\right).\tag{21.40}\label{eq:stommel_result_v} \end{align} \]

The solutions of these equations are shown for

\[ \begin{align} \beta &= 7\cdot 10^{-12}\:\frac{1}{\text{ms}},\\ D &= 4{,}000\:\text{m},\\ U &= 10\:\frac{\text{m}}{\text{s}},\\ L_x &= 10{,}000\:\text{km},\\ L_y &= 6{,}000\:\text{km},\\ \rho_0 &= 1024\:\frac{\text{kg}}{\text{m}^3},\\ \rho_a &= 1.2\:\frac{\text{kg}}{\text{m}^3},\\ C &= 0.001,\\ r &= 5.0\:\frac{\text{kg}}{\text{m}^2\text{s}} \end{align} \]

in Fig. 21.1. The Stommel model thus explains the existence of western boundary currents (Gulf Stream, Kuroshio).

In the case of the $f-$plane, $\lambda_2 = -\lambda_1$ holds. Inserting this into Eq. (21.38) yields

\[ \begin{align} \psi\left(x,y\right) &= \frac{L_y}{r\pi}C\rho_aU^2\left[\frac{1 - \exp\left(-\lambda_1L_x\right)}{\gamma}\exp\left(\lambda_1x\right) + \frac{\exp\left(\lambda_1L_x\right) - 1}{\gamma}\exp\left(-\lambda_1x\right) + 1\right]\sin\left(\pi\frac{y}{L_y}\right)\nonumber\\ \Rightarrow\psi_f\left(\frac{L_x}{2}+x,y\right) &= \frac{L_y}{r\pi}C\rho_aU^2\left[\frac{\exp\left(\lambda_1\frac{L_x}{2}\right) - \exp\left(-\lambda_1\frac{L_x}{2}\right)}{\gamma}\left(\exp\left(\lambda_1x\right) + \exp\left(-\lambda_1x\right)\right) + 1\right]\sin\left(\pi\frac{y}{L_y}\right). \end{align} \]

In this case, the solution is therefore symmetric about $x = \frac{L_x}{2}$, so no western boundary current arises. The latitude dependence of the Coriolis parameter is thus constitutive for the existence of the boundary currents.

21.3.1.2 Munk model

21.4 Thermohaline circulation

Thermohaline circulation describes the planetary circulation that arises because particles, as long as they are not subject to diabatic forcing, move along isopycnals. The basic idea is as follows:

1. Through diabatic fluxes, in particular evaporation and the resulting cooling and increase in salinity, the potential density of a water mass near the surface decreases.

2. The water mass sinks (deep-water formation), and the particles move along isopycnals through the ocean until the isopycnal reaches the surface again.

In practice, a large part of deep-water formation takes place in the North Atlantic after water transported there by the Gulf Stream has cooled through latent and sensible heat fluxes. This circulation is also known as the ocean conveyor belt. In the North Pacific, by contrast, no deep-water formation occurs because

• the trade winds transport freshwater in the form of precipitation from the Atlantic to the Pacific, so that the Atlantic is saltier, and

• water from the very salty Mediterranean reinforces deep-water formation in the Atlantic, which does not happen in the Pacific.

The water formed in the North Atlantic is called North Atlantic deep water (NADW). The densest and thus deepest water, however, is Antarctic bottom water (AABW). This water is not only very cold, but also very saline because ice formation expels salt, which gives it the highest potential density.

21.5 Ocean models

Ocean models simulate the state of the ocean, which is determined by the variables described in Sect. 21.1, analogous to atmospheric models, which simulate the state of the atmosphere. Compared to the atmosphere, ocean modeling has three major differences:

• The thermodynamic equation of state is considerably more complicated. In numerical terms, however, this is not a major problem.

• In the ocean there are no phase transitions. This is a major simplification, since it removes the need to compute phase transition rates and latent heat, and fewer tracers have to be advected.

• While the upper boundary of the atmosphere is set more or less arbitrarily, the water column has a clear vertical extent $[-H,\eta]$, where $H$ is the depth and $\eta$ is the displacement of the water surface due to currents (tides, wind-driven circulation, and waves). The difficulties here are that $H$ can become zero (land), and that $\eta$ is time dependent.

21.5.1 Prognostic variables

The prognostic variables of ocean models are essentially those described in Sect. 21.1, albeit with a few peculiarities. The following is a list of prognostic variables in typical use; other variables may of course be employed instead, as long as those listed in Sect. 21.1 can be derived from them:

• The horizontal velocities and the surface elevation are used directly as such.

• Since ocean models almost always make use of the hydrostatic approximation, the vertical velocity is expressed in the p-system or with respect to another generalized vertical coordinate.

• Instead of the temperature, the potential temperature is used, since it is a passive tracer and satisfies a simple continuity equation.

• The density is used directly as such, since, in contrast to the potential density, it satisfies a simple continuity equation.

• The salinity is used either as a specific quantity $q_s$ (mass of salt per mass of seawater) or as a volume-based quantity $\rho q_s$ (mass of salt per volume of seawater).

All other thermodynamic quantities can be derived from the potential temperature, the density, and the salinity.

21.5.2 Vertical coordinates

Due to the challenge described in the previous section, the vertical coordinate plays a particularly important role in ocean modeling. In general, a vertical coordinate $\kappa$ is a bijective mapping

\[ \begin{align} z = z\left(\kappa\right). \end{align} \]

21.5.2.1 z-coordinates

With z-coordinates (also called geopotential coordinates), the coordinate surfaces are exactly horizontal. The main advantage is the absence of metric terms in horizontal derivatives. At the coasts, the coordinate surfaces intersect the ocean floor, so that the number of layers depends on the horizontal coordinates. Owing to the surface displacement $\eta$, the thickness of the first layer is moreover time dependent. The bathymetry is thereby represented in a step-like fashion.

21.5.2.2 $\sigma_z$-coordinates

With $\sigma_z$-coordinates, the layers are fitted to the bathymetry so that they do not intersect the ocean floor. In shallow waters the layers thereby become very thin, and on slopes the coordinate surfaces run very steeply. The mathematical formulation is analogous to the terrain-following coordinates in the atmosphere described in Sect. 12.3.

21.5.2.3 $\rho$-coordinates

Within the deeper ocean, the motions proceed predominantly along isopycnals. Aligning the coordinate surfaces with these is numerically advantageous, since vertical motions then almost vanish. This is, however, only possible if the ocean is stably stratified (otherwise the potential density within a water column cannot be mapped uniquely onto the geometric height), which is not the case everywhere, which is why these coordinates are not globally applicable. This is analogous to the isentropic coordinates in the atmosphere considered in Sect. 12.2.

21.5.2.4 $z^\star$-coordinates

$z^\star$-coordinates („z-star coordinates“) are a generalization of z-coordinates that resolves the problem that the thickness of the topmost layer depends on the surface displacement. The transformation to geometric coordinates reads

\[ \begin{align} z = \eta + z^\star\frac{\eta + H}{H}. \end{align} \]

The coordinate surfaces are therefore no longer exactly horizontal, but they still intersect the ocean floor.

21.5.2.5 Hybrid coordinates

All of the vertical coordinates listed so far have substantial advantages and disadvantages. Hybrid coordinates are an attempt to combine the advantages of the various coordinates with one another. Within the upper few tens of meters, the so-called thermocline, the ocean is well mixed and weakly stratified. Here, z- or $z^\star$-coordinates are appropriate. Further down, the ocean is dominated by quasi-geostrophic motions along the isopycnals (surfaces of equal potential density). In coastal regions, by contrast, terrain-following coordinates are appropriate in order to be able to represent the bathymetry. Hybrid coordinates are accordingly a mixture of $\rho-$, $\sigma_z-$, and $z^\star-$coordinates. They give the ocean model HyCOM (Hybrid Coordinate Ocean Model) its name.

21.6 Sea state forecast

21.6.1 Prognostic variables

Sea state is another word for water surface waves. The prognostic variable considered here is the displacement of the water surface

\[ \begin{align} h = h\left(x, y, t\right) \end{align} \]

from the mean position of the water surface.$h$ cannot simply be defined as a deviation from the geoid, since the dynamic topography would still be superimposed on it. This does not contribute to the displacement associated with the waves. The averaging interval must therefore be chosen accordingly. In many situations, $h$ is not a function of the horizontal coordinates, for example in surf or in a rough sea, because in such cases the position of the water surface is no longer uniquely defined. Such effects are accounted for later in the form of energy-dissipating source terms.

One could now use the shallow-water equations derived in Sect. 13.8.1 as a set of prognostic equations, possibly supplemented by some semi-empirical additional terms. However, this has several disadvantages:

• In order to capture a significant portion of the energy of the wave field, a very small time step and a very fine spatial resolution would be required.

• In most cases, the phases of the waves are not known and are not particularly relevant. Instead, quantities such as wave direction and wave height are decisive.

Therefore, water surface waves are usually described by means of a radiative transfer equation (RTE). The prognostic variable is then the spectral radiance

\[ \begin{align} N = N\left(\mathbf{k},\mathbf{r},t\right) \end{align} \]

Here $\mathbf{r}$ is a two-dimensional position vector and $\mathbf{k}$ is a two-dimensional wave vector. Usually, $\mathbf{k}$ is given not in Cartesian coordinates but in polar coordinates $\left(k,\theta\right)$. As the dispersion relation $\omega = \omega\left(k, \theta\right)$, one uses that of water surface waves

\[ \begin{align} \omega^2 &= gk\tanh\left(kD\right) \Rightarrow \omega = \sqrt{gk\tanh\left(kD\right)}, \end{align} \]

where $D$ is the mean water depth, i.e. the water depth without waves. From this it follows that

\[ \begin{align} c_\text{ph} &= \frac{\omega}{k} = \sqrt{\frac{g\tanh\left(kD\right)}{k}},\\ c_\text{gr} &= \frac{1}{2\omega}\frac{\partial\omega^2}{\partial k} = \frac{1}{2\omega}\left[g\tanh\left(kD\right) + \frac{gkD}{\cosh^2\left(kD\right)}\right] = \frac{g\tanh\left(kD\right)}{2\omega}\left[1 + \frac{kD}{\sinh\left(kD\right)\cosh\left(kD\right)}\right]\nonumber\\ &= \frac{c_\text{ph}}{2}\left[1 + \frac{2kD}{2\sinh\left(kD\right)\cosh\left(kD\right)}\right] = \frac{c_\text{ph}}{2}\left[1 + \frac{2kD}{\sinh\left(2kD\right)}\right]. \end{align} \]

Thus the group velocity is isotropic but not homogeneous,

\[ \begin{align} c_\text{gr} = c_\text{gr}\left(k,\mathbf{r},t\right). \end{align} \]

The spatial and temporal dependence arises through the spatial and temporal dependence of $D$.

21.6.2 Wave action equation

From this, one can derive a spectral radiative flux density

\[ \begin{align} N\left(k, \theta,\mathbf{r},t\right)\left(c_\text{gr}\left(k,\mathbf{r},t\right) + \mathbf{v}\left(\mathbf{r}, t\right)\right) \end{align} \]

where $\mathbf{v} = \mathbf{v}\left(\mathbf{r}, t\right)$ is the current velocity. The radiative transfer equation is a kind of continuity equation for $N$, which is conceptually related to the conservation of energy:

The spatial and temporal dependence is no longer written explicitly here. In addition, five source terms have been included:

• $S_\text{nl}\left(k, \theta\right)$ is the energy source resulting from nonlinear effects. This term is particularly relevant for interactions across scales. One example is the emergence of long waves, the so-called swell, from the shorter-wave wind sea produced directly by the wind, via upscaling.

• $S_\text{ws}\left(k, \theta\right)$ is the energy source resulting from interaction with the wind field. This is the main energy source of the wave field.

• $S_\text{wc}\left(k, \theta\right)$ is the energy source resulting from spray, the so-called whitecapping. This is a dissipative and therefore usually negative term.

• $S_\text{diss}\left(k, \theta\right)$ is the energy source associated with internal friction, i.e. viscosity. Since viscosity acts quadratically more strongly on the smaller scales, this effect causes wind sea to be dissipated faster than swell.

• $S_\text{bd}\left(k, \theta\right)$ is the energy source associated with interaction with the bathymetry, bottom drag. This, too, is a dissipative and therefore usually negative term.

Depending on the particular situation, further source terms can be included. Eq. (21.60) is also called the wave action equation.

The total energy of the wave spectrum $E$ at a given location and time is the integral over the entire spectrum:

\[ \begin{align} E = \int_0^{2\pi}\int_0^\infty N\left(k,\theta\right)dkd\theta.\tag{21.61}\label{eq:wave_spectrum_total_energy} \end{align} \]

The energy-weighted spectral mean of a quantity $\psi$ is therefore given by

\[ \begin{align} \newoverline{\psi} \coloneqq \frac{1}{E}\int_0^{2\pi}\int_0^\infty\psi\left(k,\theta\right)N\left(k,\theta\right)dkd\theta.\tag{21.62}\label{eq:wave_spectral_average} \end{align} \]

21.6.3 Diagnostic variables

21.6.3.1 Significant wave height

The significant wave height is a kind of representative wave height for which

\[ \begin{align} H_s = 4\sqrt{E}. \end{align} \]

21.6.3.2 Mean wave direction

The mean wave direction $\theta_m$ is calculated from Eq. (21.62) as

\[ \begin{align} \theta_m = \arctan2\left(b,a\right) \end{align} \]

with

\[ \begin{align} a &\coloneqq \frac{1}{E}\int_0^{2\pi}\int_0^\infty\cos\left(\theta\right)N\left(k,\theta\right)dkd\theta,\nonumber\\ b &\coloneqq \frac{1}{E}\int_0^{2\pi}\int_0^\infty\sin\left(\theta\right)N\left(k,\theta\right)dkd\theta. \end{align} \]

21.6.3.3 Mean wavelength

For the mean wavelength $L_m$, one has

\[ \begin{align} L_m = \newoverline{\left(\frac{2\pi}{k}\right)} = 2\pi\newoverline{k^{-1}}. \end{align} \]

21.6.3.4 Mean wave period

There are different ways of calculating the mean wave period:

\[ \begin{align} T_{m,1} &= \frac{2\pi}{\newoverline{\sigma}},\\ T_{m,2} &= \frac{2\pi}{\sqrt{\newoverline{\sigma^2}}},\\ T_{m,-1} &= \newoverline{\left(\frac{2\pi}{\sigma}\right)} = 2\pi\newoverline{\sigma^{-1}}. \end{align} \]

Here $\sigma$ is the angular frequency measured relative to the seabed.

21.6.4 Sea state forecast models

The most common sea-state forecast model is Wavewatch III. Models based on Eq. (21.60) are so-called third-generation sea-state forecast models. They solve this equation on a grid and are therefore grid-point models, even though they are sometimes referred to as spectral models because their prognostic variable has spectral meaning. They solve Eq. (21.60) by introducing, at each grid point, a spectral direction-dependent grid in $\left(k,\theta\right)$ space. In most cases, a fairly large variety of semi-empirical source terms $Q_i$ is also included.