The aim of this chapter is to discretize the momentum diffusion and dissipation on the three-dimensional deformed C-grid.
First, the aim is to discretize the momentum diffusion within the framework of the shallow water equations on the two-dimensional regular hexagonal grid. In the two-dimensional case, according to Eq. (8.50), one can write for the frictional acceleration
\[ \begin{align} \frac{\partial\mathbf{v}_\text{diff}}{\partial t} = \frac{1}{h}\nabla\cdot\left[2K\left(\begin{array}{cc} S_{1, 1} & S_{1, 2}\\ S_{2, 1} & S_{2, 2} \end{array}\right)\right] = \frac{1}{h}\nabla\cdot\left[K\left(\begin{array}{cc} 2\frac{\partial v_1}{\partial x_1} & \frac{\partial v_1}{\partial x_2} + \frac{\partial v_2}{\partial x_1}\\ \frac{\partial v_2}{\partial x_1} + \frac{\partial v_1}{\partial x_2} & 2\frac{\partial v_2}{\partial x_2} \end{array}\right)\right].\tag{29.1}\label{eq:momentum_diff_hex_deriv_0} \end{align} \]
Here, the density $\rho$ was replaced by the layer thickness $h$, $K$ is a general mass-weighted diffusion coefficient. The three-dimensional divergence $\nabla\cdot\mathbf{v}$ vanishes in the incompressible case. As an ansatz for the discretization, one uses
\[ \begin{align} \frac{\partial\mathbf{v}^T_\text{diff}}{\partial t} = \left(\delta_1, \delta_2, \delta_3\right)\cdot\tau = \left(\delta_1, \delta_2, \delta_3\right)\cdot\left( \begin{array}{ccc} \tau_{1, 1} & \tau_{1, 2} & \tau_{1, 3}\\ \tau_{2, 1} & \tau_{2, 2} & \tau_{2, 3} \\ \tau_{3, 1} & \tau_{3, 2} & \tau_{3, 3} \end{array}\right) \end{align} \]
The aim of this section is to determine the $\tau_{i, j}$ uniquely. Motivated by Eq. (29.1), one requires the symmetry of $\tau$, i.e.,
\[ \begin{align} \tau_{i, j} = \tau_{j, i}. \end{align} \]
The tensor $\tau$ is developed in several modifying iterations. One first sets
\[ \begin{align} \tau = K\left( \begin{array}{ccc} 2\delta_1u_1 & \delta_2u_1 + \delta_1u_2 & \delta_3u_1 + \delta_1u_3\\ \delta_2u_1 + \delta_1u_2 & 2\delta_2u_2 & \delta_3u_2 + \delta_2u_3 \\ \delta_3u_1 + \delta_1u_3 & \delta_3u_2 + \delta_2u_3 & 2\delta_3u_3 \end{array}\right). \end{align} \]
From this it follows
\[ \begin{align} \frac{\partial u_{1, \text{diff}}}{\partial t} &= K\left[\delta_1\left(2\delta_1u_1\right) + \delta_2\left(\delta_2u_1 + \delta_1u_2\right) + \delta_3\left(\delta_3u_1 + \delta_1u_3\right)\right]\nonumber\\ &= K\left(\delta_{1, 1}u_1 + \delta_{1, 1}u_1 + \delta_{2, 2}u_1 + \delta_{3, 3}u_1 + \delta_{2, 1}u_2 + \delta_{3, 1}u_3\right)\nonumber\\ &= K\left[\delta_{1, 1}u_1 + \delta_{2, 2}u_1 + \delta_{3, 3}u_1 + \delta_1\left(\delta_1u_1 + \delta_2u_2 + \delta_3u_3\right)\right] = K\left(\frac{3}{2}\Delta u_1 + \frac{3}{2}\delta_1D\right). \end{align} \]
It follows from this
\[ \begin{align} \tau = K\left( \begin{array}{ccc} \frac{4}{3}\delta_1u_1 & \frac{2}{3}\left(\delta_2u_1 + \delta_1u_2\right) & \frac{2}{3}\left(\delta_3u_1 + \delta_1u_3\right)\\ \frac{2}{3}\left(\delta_2u_1 + \delta_1u_2\right) & \frac{4}{3}\delta_2u_2 & \frac{2}{3}\left(\delta_3u_2 + \delta_2u_3\right) \\ \frac{2}{3}\left(\delta_3u_1 + \delta_1u_3\right) & \frac{2}{3}\left(\delta_3u_2 + \delta_2u_3\right) & \frac{4}{3}\delta_3u_3 \end{array}\right) \end{align} \]
it follows for the frictional acceleration
\[ \begin{align} \frac{\partial u_{1, \text{diff}}}{\partial t} &= K\left(\Delta u_1 + \delta_1D\right). \end{align} \]
The term proportional to $\delta_1D$ does not occur in the shallow water equations, since these describe an incompressible medium, and should therefore be eliminated or at least reduced if possible. The terms that are not on the diagonal must not be modified in order not to destroy the symmetry of the tensor. Therefore, with a parameter $0 \leq \alpha \leq 1$, one sets
\[ \begin{align} \tau &= K\left( \begin{array}{ccc} \frac{4}{3}\delta_1u_1 - \alpha D & \frac{2}{3}\left(\delta_2u_1 + \delta_1u_2\right) & \frac{2}{3}\left(\delta_3u_1 + \delta_1u_3\right)\\ \frac{2}{3}\left(\delta_2u_1 + \delta_1u_2\right) & \frac{4}{3}\delta_2u_2 - \alpha D & \frac{2}{3}\left(\delta_3u_2 + \delta_2u_3\right) \\ \frac{2}{3}\left(\delta_3u_1 + \delta_1u_3\right) & \frac{2}{3}\left(\delta_3u_2 + \delta_2u_3\right) & \frac{4}{3}\delta_3u_3 - \alpha D \end{array}\right).\nonumber\\ &= K\left( \begin{array}{ccc} \frac{2}{3}\left(2 - \alpha\right)\delta_1u_1 - \alpha\frac{2}{3}\left(\delta_2u_2 + \delta_3u_3\right) & \frac{2}{3}\left(\delta_2u_1 + \delta_1u_2\right) & \frac{2}{3}\left(\delta_3u_1 + \delta_1u_3\right)\\ \frac{2}{3}\left(\delta_2u_1 + \delta_1u_2\right) & \frac{2}{3}\left(2 - \alpha\right)\delta_2u_2 - \alpha\frac{2}{3}\left(\delta_1u_1 + \delta_3u_3\right) & \frac{2}{3}\left(\delta_3u_2 + \delta_2u_3\right) \\ \frac{2}{3}\left(\delta_3u_1 + \delta_1u_3\right) & \frac{2}{3}\left(\delta_3u_2 + \delta_2u_3\right) & \frac{2}{3}\left(2 - \alpha\right)\delta_3u_3 - \alpha\frac{2}{3}\left(\delta_1u_1 + \delta_2u_2\right) \end{array}\right). \end{align} \]
From this it now follows
\[ \begin{align} \frac{\partial u_{1, \text{diff}}}{\partial t} &= K\left[\Delta u_1 + \left(1 - \alpha\right)\delta_1D\right]. \end{align} \]
The factor $1 - \alpha$ should be as close to zero as possible, so $\alpha$ should be as close to one as possible. According to Eq. (8.66), for the dissipation rate one has
\[ \begin{align} \frac{3}{2}h\epsilon &= \frac{2}{3}\left(2 - \alpha\right)\left(\delta_1u_1\right)^2 - \alpha\frac{2}{3}\left(\delta_1u_1\delta_2u_2 + \delta_1u_1\delta_3u_3\right) + \frac{2}{3}\left(\delta_1u_2\delta_2u_1 + \delta_1u_2\delta_1u_2\right) + \frac{2}{3}\left(\delta_1u_3\delta_3u_1 + \delta_1u_3\delta_1u_3\right)\nonumber\\ &+ \frac{2}{3}\left(\delta_2u_1\delta_2u_1 + \delta_2u_1\delta_1u_2\right) + \frac{2}{3}\left(2 - \alpha\right)\left(\delta_2u_2\right)^2 - \alpha\frac{2}{3}\left(\delta_2u_2\delta_1u_1 + \delta_2u_2\delta_3u_3\right) + \frac{2}{3}\left(\delta_2u_3\delta_3u_2 + \delta_2u_3\delta_2u_3\right)\nonumber\\ &+ \frac{2}{3}\left(\delta_3u_1\delta_3u_1 + \delta_3u_1\delta_1u_3\right) + \frac{2}{3}\left(\delta_3u_2\delta_3u_2 + \delta_3u_2\delta_2u_3\right) + \frac{2}{3}\left(2 - \alpha\right)\left(\delta_3u_3\right)^2 - \alpha\frac{2}{3}\left(\delta_3u_3\delta_1u_1 + \delta_3u_3\delta_2u_2\right)\nonumber\\ &= \frac{2}{3}\left(2 - \alpha\right)\left(\delta_1u_1\right)^2 - \alpha\frac{2}{3}\left(\delta_1u_1\delta_2u_2 + \delta_1u_1\delta_3u_3\right) + \textcolor{red}{\frac{2}{3}\left(\delta_2u_1 + \delta_1u_2\right)^2} + \textcolor{red}{\frac{2}{3}\left(\delta_3u_1 + \delta_1u_3\right)^2}\nonumber\\ &+ \frac{2}{3}\left(2 - \alpha\right)\left(\delta_2u_2\right)^2 - \alpha\frac{2}{3}\left(\delta_2u_2\delta_1u_1 + \delta_2u_2\delta_3u_3\right) + \textcolor{red}{\frac{2}{3}\left(\delta_3u_2 + \delta_2u_3\right)^2}\nonumber\\ &+ \frac{2}{3}\left(2 - \alpha\right)\left(\delta_3u_3\right)^2 - \alpha\frac{2}{3}\left(\delta_3u_3\delta_1u_1 + \delta_3u_3\delta_2u_2\right). \end{align} \]
The terms marked in red arise from the off-diagonal elements and are all squares and therefore non-negative. For the remaining terms one has
\[ \begin{align} &\frac{2}{3}\left(2 - \alpha\right)\left(\delta_1u_1\right)^2 - \alpha\frac{2}{3}\left(\delta_1u_1\delta_2u_2 + \delta_1u_1\delta_3u_3\right) + \frac{2}{3}\left(2 - \alpha\right)\left(\delta_2u_2\right)^2 - \alpha\frac{2}{3}\left(\delta_2u_2\delta_1u_1 + \delta_2u_2\delta_3u_3\right)\nonumber\\ &+ \frac{2}{3}\left(2 - \alpha\right)\left(\delta_3u_3\right)^2 - \alpha\frac{2}{3}\left(\delta_3u_3\delta_1u_1 + \delta_3u_3\delta_2u_2\right)\nonumber\\ &= \frac{2}{3}\left(2 - \alpha\right)\left(\delta_1u_1\right)^2 + \frac{2}{3}\left(2 - \alpha\right)\left(\delta_2u_2\right)^2 + \frac{2}{3}\left(2 - \alpha\right)\left(\delta_3u_3\right)^2 - \alpha\frac{4}{3}\left(\delta_1u_1\delta_2u_2 + \delta_1u_1\delta_3u_3 + \delta_2u_2\delta_3u_3\right)\nonumber\\ &\hastobe \alpha\frac{2}{3}\left(\delta_1u_1 - \delta_2u_2\right)^2 + \alpha\frac{2}{3}\left(\delta_1u_1 - \delta_3u_3\right)^2 + \alpha\frac{2}{3}\left(\delta_2u_2 - \delta_3u_3\right)^2. \end{align} \]
The requirement in the last step is motivated by the wish to obtain a non-negative dissipation rate. This yields
\[ \begin{align} \frac{2}{3}\left(2 - \alpha\right) = \frac{4}{3}\alpha \Leftrightarrow \frac{4}{3} - \alpha\frac{2}{3} = \frac{4}{3}\alpha \Rightarrow \frac{4}{3} = \frac{6}{3}\alpha \Leftrightarrow\alpha = \frac{2}{3}. \end{align} \]
With this choice of $\alpha$, the terms on the main diagonal also contribute non-negatively to the dissipation rate and one has
\[ \begin{align} \epsilon \geq 0. \end{align} \]
The result for the friction tensor thus reads
\[ \begin{align} \tau &= K\left( \begin{array}{ccc} \frac{4}{9}\left[2\delta_1u_1 - \left(\delta_2u_2 + \delta_3u_3\right)\right] & \frac{2}{3}\left(\delta_2u_1 + \delta_1u_2\right) & \frac{2}{3}\left(\delta_3u_1 + \delta_1u_3\right)\\ \frac{2}{3}\left(\delta_2u_1 + \delta_1u_2\right) & \frac{4}{9}\left[2\delta_2u_2 - \left(\delta_1u_1 + \delta_3u_3\right)\right] & \frac{2}{3}\left(\delta_3u_2 + \delta_2u_3\right)\\ \frac{2}{3}\left(\delta_3u_1 + \delta_1u_3\right) & \frac{2}{3}\left(\delta_3u_2 + \delta_2u_3\right) & \frac{4}{9}\left[2\delta_3u_3 - \left(\delta_1u_1 + \delta_2u_2\right)\right] \end{array}\right). \end{align} \]
The generalization
\[ \begin{align} \tau &= \left( \begin{array}{ccc} K_c\frac{4}{9}\left[2\delta_1u_1 - \left(\delta_2u_2 + \delta_3u_3\right)\right] & K_1\frac{2}{3}\left(\delta_2u_1 + \delta_1u_2\right) & K_2\frac{2}{3}\left(\delta_3u_1 + \delta_1u_3\right)\\ K_1\frac{2}{3}\left(\delta_2u_1 + \delta_1u_2\right) & K_c\frac{4}{9}\left[2\delta_2u_2 - \left(\delta_1u_1 + \delta_3u_3\right)\right] & K_3\frac{2}{3}\left(\delta_3u_2 + \delta_2u_3\right)\\ K_2\frac{2}{3}\left(\delta_3u_1 + \delta_1u_3\right) & K_3\frac{2}{3}\left(\delta_3u_2 + \delta_2u_3\right) & K_c\frac{4}{9}\left[2\delta_3u_3 - \left(\delta_1u_1 + \delta_2u_2\right)\right] \end{array}\right)\tag{29.15}\label{eq:stress_tensor_hex_regular} \end{align} \]
does not destroy any of the verified statements.
The divergence of the tensor Eq. (29.15) does not yet satisfy the Thuburn condition. Motivated by Sect. 27.6.1, one therefore makes the following ansatz for the frictional acceleration at the edge $e$:
\[ \begin{align} h_e\frac{\partial u_{e, \text{diff}}}{\partial t} &= \left[\nabla\left(K_cE_{c(e)}\right)\right]_e + \frac{1}{l_e}\left(\frac{1}{3}\sum_{r\in t_1}K_rF_{r(e)} - \frac{1}{3}\sum_{r\in t_2}K_rF_{r(e)}\right)\nonumber\\ &= \left[\nabla\left(K_cE_{c(e)}\right)\right]_e + \frac{1}{3l_e}\left(\sum_{r\in t_1}K_rF_{r(e)} - \sum_{r\in t_2}K_rF_{r(e)}\right). \end{align} \]
Here $E_c$ and $F_r$ are deformations, which are localized in the cell centers and at the edges, respectively. The $E_c$ are located on the diagonal of the stress tensor; they are referred to as strain deformations. The $F_r$ are not located on the diagonal of the stress tensor; they are referred to as shear deformations. Here one allows both to depend not only on their position but also on the edge at which one wants to compute the acceleration, i.e. $E_{c(e)}$ and $F_{r(e)}$, respectively. The triangle $t_1$ is located with respect to the edge $e$ in the direction $\mathbf{k}\times\mathbf{n}_e$, the triangle $t_2$ in the direction $-\mathbf{k}\times\mathbf{n}_e$. $K_c$ and $K_r$ are non-negative diffusion coefficients, which are left open in this chapter. Without loss of generality, one assumes that $\mathbf{n}_e$ is aligned parallel to $\mathbf{i}_1$.
On the regular grid one obtains for the strain deformation
\[ \begin{align} E_{c(e)} &= \frac{4}{9}\left[2\delta_1u_1 - \left(\delta_2u_2 + \delta_3u_3\right)\right] = \frac{4}{9}\left[3\delta_1u_1 - \left(\delta_1u_1 + \delta_2u_2 + \delta_3u_3\right)\right] = \frac{4}{3}\delta_1u_1 - \frac{4}{9}\left(\delta_1u_1 + \delta_2u_2 + \delta_3u_3\right)\nonumber\\ &= -\frac{4}{9}\left(\delta_1u_1 + \delta_2u_2 + \delta_3u_3\right) + \frac{4}{3}\delta_1u_1 = -\frac{2}{3}\frac{2}{3}\left(\delta_1u_1 + \delta_2u_2 + \delta_3u_3\right) + \frac{4}{3}\delta_1u_1\nonumber\\ &\stackrel{\href{ch-26-problems-of-a-three-element-generating-s.html#eq:div_3elements2d}{\text{Eq. (27.33)}}}{=} -\frac{2}{3}D + \frac{4}{3}\delta_1u_1 \end{align} \]
For the shear deformations, according to Eq. (29.15)
\[ \begin{align} &\delta_2\left[\frac{2}{3}\left(\delta_2u_1 + \delta_1u_2\right)\right] + \delta_3\left[\frac{2}{3}\left(\delta_3u_1 + \delta_1u_3\right)\right] \hastobe \frac{1}{3}\delta_y\left(F_{1(e)} + F_{2(e)} + F_{3(e)}\right)\nonumber\\ &\stackrel{\href{ch-26-problems-of-a-three-element-generating-s.html#eq:ddy_hex}{\text{Eq. (27.31)}}}{=} \frac{1}{3\sqrt{3}}\left(\delta_2 - \delta_3\right)\left(F_{1(e)} + F_{2(e)} + F_{3(e)}\right)\tag{29.18}\label{eq:deriv_hex_diff_deformed_0} \end{align} \]
must hold. The choice of $F_{1(e)}$ is arbitrary here, since this term is first added and then subtracted again. One therefore naturally chooses
\[ \begin{align} F_{1(e)} &= \zeta_1. \end{align} \]
First, one treats Eq. (29.18) in the form
\[ \begin{align} &\delta_2\left[\frac{2}{3}\left(\delta_2u_1 + \delta_1u_2\right)\right] + \delta_3\left[\frac{2}{3}\left(\delta_3u_1 + \delta_1u_3\right)\right] = \frac{1}{\sqrt{3}}\left(\delta_2F_{3(e)} - \delta_3F_{2(e)}\right). \end{align} \]
With this, one obtains
\[ \begin{align} F_{2(e)} &= -\frac{2}{\sqrt{3}}\left(\delta_3u_1 + \delta_1u_3\right) = \frac{2}{\sqrt{3}}\left(-\delta_3u_1 - \delta_1u_3\right) = \frac{2}{\sqrt{3}}\left[\left(\delta_3u_1 - \delta_1u_3\right) - 2\delta_3u_1\right] \stackrel{\href{ch-26-problems-of-a-three-element-generating-s.html#eq:hex_curl_rhombus_1}{\text{Eq. (27.171)}}}{=} \zeta_2 - \frac{4}{\sqrt{3}}\delta_3u_1,\\ F_{3(e)} &= \frac{2}{\sqrt{3}}\left(\delta_2u_1 + \delta_1u_2\right) = \frac{2}{\sqrt{3}}\left[\left(\delta_1u_2 - \delta_2u_1\right) + 2\delta_2u_1\right] \stackrel{\href{ch-26-problems-of-a-three-element-generating-s.html#eq:hex_curl_rhombus_2}{\text{Eq. (27.172)}}}{=} \zeta_3 + \frac{4}{\sqrt{3}}\delta_2u_1. \end{align} \]
Analogously to Sect. 27.6.1, one now uses the averaged shear deformations, which leads to Eq. (29.18). The procedure for calculating the gradient of the component $u_1$ is presented in Sect. 29.2.0.1 in the context of the deformed grid.
In this section the results of the previous section are generalized to the deformed hexagonal C-grid. Here too one remains within the framework of the shallow water equations. Essentially, one can use the same formulas as on the regular grid, but the velocity gradients can be written down slightly differently. Because of
\[ \begin{align} F_{2(e)} + F_{3(e)} &= \zeta_2 - \frac{4}{\sqrt{3}}\delta_3u_1 + \zeta_3 + \frac{4}{\sqrt{3}}\delta_2u_1 = \zeta_2 + \zeta_3 + 4\frac{1}{\sqrt{3}}\left(\delta_2 - \delta_3\right)u_1 = \zeta_2 + \zeta_3 + 4\delta_yu_1\nonumber\\ &= \zeta_2 + 2\delta_yu_1 + \zeta_3 + 2\delta_yu_1 \end{align} \]
one can also write the deformations in the form
\[ \begin{align} E_{c(e)} &= -\frac{2}{3}D + \frac{4}{3}\delta_xu_1\vert_c,\\ F_{1(e)} &= \zeta_1,\\ F_{2(e)} &= \zeta_2 + 2\delta_yu_1\vert_2,\\ F_{3(e)} &= \zeta_3 + 2\delta_yu_1\vert_3 \end{align} \]
The velocity gradients are calculated over cells and rhombuses, respectively.
For the dissipation rate $\epsilon_\text{hex}$, one can simply set
\[ \begin{align} \epsilon_\text{hex} = -\rho\frac{\partial\mathbf{v}_\text{diff}}{\partial t}\cdot\mathbf{v} \end{align} \]
Although this does not correspond exactly to the molecular dissipation rate Eq. (8.66), the global integrals over both ansätze are identical according to Eq. (8.96). Thus, energetic self-consistency is ensured.