The two unit vectors of the plane are denoted by
\[ \begin{align} \mathbf{i} &\coloneqq \left(\begin{array}{c} 1\\ 0 \end{array}\right),\\ \mathbf{j} &\coloneqq \left(\begin{array}{c} 0\\ 1 \end{array}\right)\\ \end{align} \]
Now one further defines a three-element generating system $\left(\mathbf{i}_1, \mathbf{i}_2 \mathbf{i}_3\right)$ by
\[ \begin{align} \mathbf{i}_1 &\coloneqq \left(\begin{array}{c} 1\\ 0 \end{array}\right) = \mathbf{i},\\ \mathbf{i}_2 &\coloneqq \left(\begin{array}{c} -\sin\left(30^\circ\right)\\ \cos\left(30^\circ\right) \end{array}\right) = \left(\begin{array}{c} -\frac{1}{2}\\ \frac{\sqrt{3}}{2} \end{array}\right) = -\frac{1}{2}\mathbf{i} + \frac{\sqrt{3}}{2}\mathbf{j},\\ \mathbf{i}_3 &\coloneqq \left(\begin{array}{c} -\frac{1}{2}\\ -\frac{\sqrt{3}}{2} \end{array}\right) = -\frac{1}{2}\mathbf{i} - \frac{\sqrt{3}}{2}\mathbf{j}, \end{align} \]
whose elements are each rotated by 120$^\circ$ relative to one another. We further define a generating system $\left(\mathbf{j}_1, \mathbf{j}_2 \mathbf{j}_3\right)$ rotated by 90$^\circ$ relative to it:
\[ \begin{align} \mathbf{j}_1 &\coloneqq \mathbf{k}\times\mathbf{i}_1 = \left(\begin{array}{c} 0\\ 0\\ 1 \end{array}\right)\times\left(\begin{array}{c} 1\\ 0\\ 0 \end{array}\right) = \left(\begin{array}{c} 0\\ 1 \end{array}\right),\\ \mathbf{j}_2 &\coloneqq \mathbf{k}\times\mathbf{i}_2 = \left(\begin{array}{c} 0\\ 0\\ 1 \end{array}\right)\times\left(\begin{array}{c} -\frac{1}{2}\\ \frac{\sqrt{3}}{2}\\ 0 \end{array}\right) = \left(\begin{array}{c} -\frac{\sqrt{3}}{2}\\ -\frac{1}{2} \end{array}\right) = -\frac{\sqrt{3}}{2}\mathbf{i} - \frac{1}{2}\mathbf{j},\\ \mathbf{j}_3 &\coloneqq \mathbf{k}\times\mathbf{i}_3 = \left(\begin{array}{c} 0\\ 0\\ 1 \end{array}\right)\times\left(\begin{array}{c} -\frac{1}{2}\\ -\frac{\sqrt{3}}{2}\\ 0 \end{array}\right) = \left(\begin{array}{c} \frac{\sqrt{3}}{2}\\ -\frac{1}{2} \end{array}\right) = \frac{\sqrt{3}}{2}\mathbf{i} - \frac{1}{2}\mathbf{j}. \end{align} \]
One observes
\[ \begin{align} \mathbf{i}_1 + \mathbf{i}_2 + \mathbf{i}_3 = \mathbf{0}, & {} & \mathbf{j}_1 + \mathbf{j}_2 + \mathbf{j}_3 = \mathbf{0}. \end{align} \]
For a two-dimensional vector $\mathbf{v}$ one can write
\[ \begin{align} \mathbf{v} = u\mathbf{i} + v\mathbf{j} \end{align} \]
where
\[ \begin{align} u = \mathbf{i}\cdot\mathbf{v}, & {} & v = \mathbf{j}\cdot\mathbf{v}. \end{align} \]
Since the $\mathbf{i}_k$, $\mathbf{j}_l$ are each pairwise linearly independent, one can write
\[ \begin{align} \mathbf{v} = u_k'\mathbf{i}_k + u_l'\mathbf{i}_l = v_k'\mathbf{j}_k + v_l'\mathbf{j}_l. \end{align} \]
Since the choice of $k$, $l$ is not unique and one could also use all three unit vectors, one can introduce another linear condition. One uses this freedom to require the notation
\[ \begin{align} \mathbf{v} &= \frac{2}{3}\left(u_1\mathbf{i}_1 + u_2\mathbf{i}_2 + u_3\mathbf{i}_3\right) = \frac{2}{3}\left(v_1\mathbf{j}_1 + v_2\mathbf{j}_2 + v_3\mathbf{j}_3\right) \end{align} \]
This is satisfied by
\[ \begin{align} u_k = \mathbf{i}_k\cdot\mathbf{v}, & {} & v_k = \mathbf{j}_k\cdot\mathbf{v} \end{align} \]
since it follows from this
\[ \begin{align} \mathbf{v} &= u\mathbf{i} + v\mathbf{j} = \frac{2}{3}\left(\frac{6}{4}u\mathbf{i} + \frac{6}{4}v\mathbf{j}\right)\nonumber\\ &= \frac{2}{3}\left(u\mathbf{i} + \frac{1}{4}u\mathbf{i} - \frac{\sqrt{3}}{4}u\mathbf{j} - \frac{\sqrt{3}}{4}v\mathbf{i} + \frac{3}{4}v\mathbf{j} + \frac{1}{4}u\mathbf{i} + \frac{\sqrt{3}}{4}u\mathbf{j} + \frac{\sqrt{3}}{4}v\mathbf{i} + \frac{3}{4}v\mathbf{j}\right)\nonumber\\ &= \frac{2}{3}\left(u\mathbf{i} + \left(-\frac{1}{2}u + \frac{\sqrt{3}}{2}v\right)\left(-\frac{1}{2}\mathbf{i} + \frac{\sqrt{3}}{2}\mathbf{j}\right) + \left(-\frac{1}{2}u - \frac{\sqrt{3}}{2}v\right)\left(-\frac{1}{2}\mathbf{i} - \frac{\sqrt{3}}{2}\mathbf{j}\right)\right)\nonumber\\ &= \frac{2}{3}\left(u_1\mathbf{i}_1 + u_2\mathbf{i}_2 + u_3\mathbf{i}_3\right), \end{align} \]
where
\[ \begin{align} u_1 &= u, & {} & u_2 = -\frac{1}{2}u + \frac{\sqrt{3}}{2}v, & {} & u_3 = -\frac{1}{2}u - \frac{\sqrt{3}}{2}v. \end{align} \]
was used. This applies analogously to the $v_k$. From this it follows
and analogously for the $v_k$. For the gradient $\nabla\alpha$ of a scalar field $\alpha$ one now obtains
\[ \begin{align} \nabla\alpha = \frac{2}{3}\left[\left(\mathbf{i}_1\cdot\nabla\alpha\right)\mathbf{i}_1 + \left(\mathbf{i}_2\cdot\nabla\alpha\right)\mathbf{i}_2 + \left(\mathbf{i}_3\cdot\nabla\alpha\right)\mathbf{i}_3\right]. \end{align} \]
Because of
\[ \begin{align} \mathbf{i}_k\cdot\nabla\alpha = \frac{\partial\alpha}{\partial x_k} \end{align} \]
it follows
With this one furthermore establishes
It holds that
\[ \begin{align} \mathbf{i}_1 &= \frac{1}{\sqrt{3}}\left(\mathbf{j}_3 - \mathbf{j}_2\right) \end{align} \]
and cyclically:
\[ \begin{align} \mathbf{i}_2 = \frac{1}{\sqrt{3}}\left(\mathbf{j}_1 - \mathbf{j}_3\right), & {} & \mathbf{i}_3 = \frac{1}{\sqrt{3}}\left(\mathbf{j}_2 - \mathbf{j}_1\right) \end{align} \]
From this it follows
\[ \begin{align} \nabla\alpha = \frac{2}{3\sqrt{3}}\left(\frac{\partial\alpha}{\partial x_1}\left(\mathbf{j}_3 - \mathbf{j}_2\right) + \frac{\partial\alpha}{\partial x_2}\left(\mathbf{j}_1 - \mathbf{j}_3\right) + \frac{\partial\alpha}{\partial x_3}\left(\mathbf{j}_2 - \mathbf{j}_1\right)\right). \end{align} \]
Rearranging yields
\[ \begin{align} \nabla\alpha = \frac{2}{3\sqrt{3}}\left(\left(\frac{\partial\alpha}{\partial x_2} - \frac{\partial\alpha}{\partial x_3}\right)\mathbf{j}_1 + \left(\frac{\partial\alpha}{\partial x_3} - \frac{\partial\alpha}{\partial x_1}\right)\mathbf{j}_2 + \left(\frac{\partial\alpha}{\partial x_1} - \frac{\partial\alpha}{\partial x_2}\right)\mathbf{j}_3\right). \end{align} \]
One has
\[ \begin{align} \mathbf{v} &= \frac{2}{3}\left(u_1\mathbf{i}_1 + u_2\mathbf{i}_2 + u_3\mathbf{i}_3\right) = \frac{2}{3}\left[u_1\left(\begin{array}{c}1\\0\end{array}\right) + u_2\left(\begin{array}{c}-\frac{1}{2}\\\frac{\sqrt{3}}{2}\end{array}\right) + u_3\left(\begin{array}{c}-\frac{1}{2}\\-\frac{\sqrt{3}}{2}\end{array}\right)\right] = \frac{2}{3}\left(\begin{array}{c}u_1 - \frac{u_2}{2} - \frac{u_3}{2}\\\frac{\sqrt{3}}{2}u_2 - \frac{\sqrt{3}}{2}u_3\end{array}\right)\nonumber\\ &= \left(\begin{array}{c}\frac{2}{3}u_1 - \frac{u_2}{3} - \frac{u_3}{3}\\\frac{1}{\sqrt{3}}u_2 - \frac{1}{\sqrt{3}}u_3\end{array}\right). \end{align} \]
From this one obtains for the divergence
\[ \begin{align} D \coloneqq \nabla\cdot\mathbf{v} = \frac{2}{3}\frac{u_1}{\partial x} - \frac{1}{3}\frac{u_2}{\partial x} - \frac{1}{3}\frac{u_3}{\partial x} + \frac{1}{\sqrt{3}}\frac{u_2}{\partial y} - \frac{1}{\sqrt{3}}\frac{u_3}{\partial y}. \end{align} \]
One has
\[ \begin{align} \frac{\partial}{\partial x} = \frac{\partial}{\partial x_1}, \end{align} \]
because of
\[ \begin{align} \mathbf{j} = \frac{1}{\sqrt{3}}\left(\mathbf{i}_2 - \mathbf{i}_3\right) \end{align} \]
one furthermore has
\[ \begin{align} \frac{\partial}{\partial y} = \frac{1}{\sqrt{3}}\left(\frac{\partial}{\partial x_2} - \frac{\partial}{\partial x_3}\right).\tag{31.31}\label{eq:ddy_hex} \end{align} \]
From this it follows
\[ \begin{align} D &= \frac{2}{3}\frac{\partial u_1}{\partial x} - \frac{1}{3}\frac{\partial u_2}{\partial x} - \frac{1}{3}\frac{\partial u_3}{\partial x} + \frac{1}{3}\left(\frac{\partial u_2}{\partial x_2} - \frac{\partial u_2}{\partial x_3}\right) - \frac{1}{3}\left(\frac{\partial u_3}{\partial x_2} - \frac{\partial u_3}{\partial x_3}\right)\nonumber\\ &= \frac{2}{3}\frac{\partial u_1}{\partial x_1} + \frac{1}{3}\frac{u_2}{\partial x_2} + \frac{1}{3}\frac{u_3}{\partial x_3} - \frac{1}{3}\frac{u_2}{\partial x_1} - \frac{1}{3}\frac{\partial u_3}{\partial x_1} - \frac{1}{3}\frac{\partial u_2}{\partial x_3} - \frac{1}{3}\frac{\partial u_3}{\partial x_2}\nonumber\\ &= \frac{2}{3}\left(\frac{\partial u_1}{\partial x_1} + \frac{\partial u_2}{\partial x_2} + \frac{\partial u_3}{\partial x_3}\right) - \frac{1}{3}\left(\frac{\partial u_2}{\partial x_1} + \frac{\partial u_2}{\partial x_2} + \frac{\partial u_2}{\partial x_3} + \frac{\partial u_3}{\partial x_1} + \frac{\partial u_3}{\partial x_2} + \frac{\partial u_3}{\partial x_3}\right). \end{align} \]
With Eq. (31.22) follows
Combining with Eq. (31.21), one further obtains
\[ \begin{align} \Delta\alpha = \frac{2}{3}\left(\frac{\partial^2\alpha}{\partial x_1^2} + \frac{\partial^2\alpha}{\partial x_2^2} + \frac{\partial^2\alpha}{\partial x_3^2}\right). \end{align} \]
One has
\[ \begin{align} u_1 &= \mathbf{v}\cdot\mathbf{i}_1 = \frac{2}{3}\left(v_1\mathbf{j}_1\cdot\mathbf{i}_1 + v_2\mathbf{j}_2\cdot\mathbf{i}_1 + v_3\mathbf{j}_3\cdot\mathbf{i}_1\right) = \frac{2}{3}\left(0\cdot v_1 - \frac{\sqrt{3}}{2}v_2 + \frac{\sqrt{3}}{2}v_3\right)\nonumber\\ &= \frac{1}{\sqrt{3}}\left(v_3 - v_2\right). \end{align} \]
By cyclically shifting the indices one obtains
\[ \begin{align} u_2 = \frac{1}{\sqrt{3}}\left(v_1 - v_3\right), & {} & u_3 = \frac{1}{\sqrt{3}}\left(v_2 - v_1\right). \end{align} \]
This gives another representation of the divergence:
\[ \begin{align} D = \frac{2}{3\sqrt{3}}\left[\left(\frac{\partial v_2}{\partial x_3} - \frac{\partial v_3}{\partial x_2}\right) + \left(\frac{\partial v_3}{\partial x_1} - \frac{\partial v_1}{\partial x_3}\right) + \left(\frac{\partial v_1}{\partial x_2} - \frac{\partial v_2}{\partial x_1}\right)\right]\tag{31.37}\label{eq:div_3elements2d_mod} \end{align} \]
One has
\[ \begin{align} \frac{\partial v_2}{\partial x_3} - \frac{\partial v_3}{\partial x_2} &= \frac{\partial v_2}{\partial x_3} + \frac{\partial v_2}{\partial x_2} - \frac{\partial v_2}{\partial x_2} - \frac{\partial v_3}{\partial x_2} = \left(\frac{\partial}{\partial x_3} + \frac{\partial}{\partial x_2}\right)v_2 - \frac{\partial}{\partial x_2}\left(v_2 + v_3\right)\nonumber\\ &= -\frac{\partial v_2}{\partial x_1} + \frac{\partial v_1}{\partial x_2} = \frac{\partial v_1}{\partial x_2} - \frac{\partial v_2}{\partial x_1}. \end{align} \]
By cyclically permuting this equation, one sees that all three bracketed summands in Eq. (31.37) are equal. From this follow
\[ \begin{align} D &= \frac{2}{\sqrt{3}}\left(\frac{\partial v_2}{\partial x_3} - \frac{\partial v_3}{\partial x_2}\right),\\ D &= \frac{2}{\sqrt{3}}\left(\frac{\partial v_3}{\partial x_1} - \frac{\partial v_1}{\partial x_3}\right),\\ D &= \frac{2}{\sqrt{3}}\left(\frac{\partial v_1}{\partial x_2} - \frac{\partial v_2}{\partial x_1}\right). \end{align} \]
One has
\[ \begin{align} \mathbf{v} &= \frac{2}{3}\left(v_1\mathbf{j}_1 + v_2\mathbf{j}_2 + v_3\mathbf{j}_3\right) = \frac{2}{3}\left[v_1\left(\begin{array}{c}0\\1\end{array}\right) + v_2\left(\begin{array}{c}-\frac{\sqrt{3}}{2}\\-\frac{1}{2}\end{array}\right) + v_3\left(\begin{array}{c}\frac{\sqrt{3}}{2}\\-\frac{1}{2}\end{array}\right)\right]\nonumber\\ &= \left(\begin{array}{c}\frac{1}{\sqrt{3}}v_3 - \frac{1}{\sqrt{3}}v_2\\\frac{2}{3}v_1 - \frac{v_2}{3} - \frac{v_3}{3}\end{array}\right). \end{align} \]
For the vorticity one obtains
\[ \begin{align} \zeta &= \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y} = \frac{\partial}{\partial x}\left(\frac{2}{3}v_1 - \frac{v_2}{3} - \frac{v_3}{3}\right) - \frac{\partial}{\partial y}\left(\frac{1}{\sqrt{3}}v_3 - \frac{1}{\sqrt{3}}v_2\right)\nonumber\\ &= \frac{\partial}{\partial x_1}\left(\frac{2}{3}v_1 - \frac{v_2}{3} - \frac{v_3}{3}\right) - \frac{1}{\sqrt{3}}\left(\frac{\partial}{\partial x_2} - \frac{\partial}{\partial x_3}\right)\left(\frac{1}{\sqrt{3}}v_3 - \frac{1}{\sqrt{3}}v_2\right)\nonumber\\ &= \frac{2}{3}\frac{\partial v_1}{\partial x_1} + \frac{1}{3}\frac{\partial v_2}{\partial x_2} + \frac{1}{3}\frac{\partial v_3}{\partial x_3} - \frac{1}{3}\frac{\partial v_2}{\partial x_1} - \frac{1}{3}\frac{\partial v_3}{\partial x_1} - \frac{1}{3}\frac{\partial v_3}{\partial x_2} - \frac{1}{3}\frac{\partial v_2}{\partial x_3}\nonumber\\ &= \frac{2}{3}\left(\frac{\partial v_1}{\partial x_1} + \frac{\partial v_2}{\partial x_2} + \frac{\partial v_3}{\partial x_3}\right) - \frac{1}{3}\left(\frac{\partial v_2}{\partial x_1} + \frac{\partial v_2}{\partial x_2} + \frac{\partial v_2}{\partial x_3} + \frac{\partial v_3}{\partial x_1} + \frac{\partial v_3}{\partial x_2} + \frac{\partial v_3}{\partial x_3}\right). \end{align} \]
With Eq. (31.22) follows
\[ \begin{align} \zeta = \frac{2}{3}\left(\frac{\partial v_1}{\partial x_1} + \frac{\partial v_2}{\partial x_2} + \frac{\partial v_3}{\partial x_3}\right). \end{align} \]
One has
\[ \begin{align} v_1 &= \mathbf{v}\cdot\mathbf{j}_1 = \frac{2}{3}\left(u_1\mathbf{i}_1 + u_2\mathbf{i}_2 + u_3\mathbf{i}_3\right)\cdot\mathbf{j}_1 = \frac{2}{3}\left(u_1\mathbf{i}_1 + u_2\mathbf{i}_2 + u_3\mathbf{i}_3\right)\cdot\mathbf{j}\nonumber\\ &= \frac{2}{3}\left(\frac{\sqrt{3}}{2}u_2 - \frac{\sqrt{3}}{2}u_3\right) = \frac{1}{\sqrt{3}}\left(u_2 - u_3\right).\tag{31.45}\label{eq:turn_3elemets2d} \end{align} \]
By cyclically shifting the indices one obtains
\[ \begin{align} v_2 = \frac{1}{\sqrt{3}}\left(u_3 - u_1\right), & {} & v_3 = \frac{1}{\sqrt{3}}\left(u_1 - u_2\right). \end{align} \]
From this follows another representation of the vorticity:
\[ \begin{align} \zeta = \frac{2}{3\sqrt{3}}\left(\frac{\partial u_2}{\partial x_1} - \frac{\partial u_3}{\partial x_1} + \frac{\partial u_3}{\partial x_2} - \frac{\partial u_1}{\partial x_2} + \frac{\partial u_1}{\partial x_3} - \frac{\partial u_2}{\partial x_3}\right) \end{align} \]
By rearranging one obtains
\[ \begin{align} \zeta = \frac{2}{3\sqrt{3}}\left[\left(\frac{\partial u_3}{\partial x_2} - \frac{\partial u_2}{\partial x_3}\right) + \left(\frac{\partial u_1}{\partial x_3} - \frac{\partial u_3}{\partial x_1}\right) + \left(\frac{\partial u_2}{\partial x_1} - \frac{\partial u_1}{\partial x_2}\right)\right]. \end{align} \]
Analogously to the case of Eq. (31.37), the following hold here as well
\[ \begin{align} \zeta &= \frac{2}{\sqrt{3}}\left(\frac{\partial u_3}{\partial x_2} - \frac{\partial u_2}{\partial x_3}\right),\tag{31.49}\label{eq:hex_curl_symm_0}\\ \zeta &= \frac{2}{\sqrt{3}}\left(\frac{\partial u_1}{\partial x_3} - \frac{\partial u_3}{\partial x_1}\right),\tag{31.50}\label{eq:hex_curl_symm_1}\\ \zeta &= \frac{2}{\sqrt{3}}\left(\frac{\partial u_2}{\partial x_1} - \frac{\partial u_1}{\partial x_2}\right).\tag{31.51}\label{eq:hex_curl_symm_2} \end{align} \]
Using the fundamental theorem of vector analysis one can write the vector field $\mathbf{v}$ in the form
\[ \begin{align} \mathbf{v} = \mathbf{k}\times\nabla\psi + \nabla\chi \end{align} \]
where
\[ \begin{align} \mathbf{k}\times\nabla\psi &= \frac{2}{3}\left(\frac{\partial\psi}{\partial x_1}\mathbf{j}_1 + \frac{\partial\psi}{\partial x_2}\mathbf{j}_2 + \frac{\partial\psi}{\partial x_3}\mathbf{j}_3\right),\\ \nabla\chi &= \frac{2}{3}\left(\frac{\partial\chi}{\partial x_1}\mathbf{i}_1 + \frac{\partial\chi}{\partial x_2}\mathbf{i}_2 + \frac{\partial\chi}{\partial x_3}\mathbf{i}_3\right). \end{align} \]
From this follow
\[ \begin{align} u_1 &= \frac{1}{\sqrt{3}}\left(\frac{\partial\psi}{\partial x_3} - \frac{\partial\psi}{\partial x_2}\right) + \frac{\partial\chi}{\partial x_1} = -\frac{1}{\sqrt{3}}\left(\frac{\partial\psi}{\partial x_2} - \frac{\partial\psi}{\partial x_3}\right) + \frac{\partial\chi}{\partial x_1}\tag{31.55}\label{eq:trivariate_helmholtz_eq_0},\\ u_2 &= -\frac{1}{\sqrt{3}}\left(\frac{\partial\psi}{\partial x_3} - \frac{\partial\psi}{\partial x_1}\right) + \frac{\partial\chi}{\partial x_2}\tag{31.56}\label{eq:trivariate_helmholtz_eq_1},\\ u_3 &= -\frac{1}{\sqrt{3}}\left(\frac{\partial\psi}{\partial x_1} - \frac{\partial\psi}{\partial x_2}\right) + \frac{\partial\chi}{\partial x_3}\tag{31.57}\label{eq:trivariate_helmholtz_eq_2}. \end{align} \]
and
\[ \begin{align} v_1 &= \frac{\partial\psi}{\partial x_1} + \frac{1}{\sqrt{3}}\left(\frac{\partial\chi}{\partial x_2} - \frac{\partial\chi}{\partial x_3}\right)\tag{31.58}\label{eq:trivariate_helmholtz_eq_3},\\ v_2 &= \frac{\partial\psi}{\partial x_2} + \frac{1}{\sqrt{3}}\left(\frac{\partial\chi}{\partial x_3} - \frac{\partial\chi}{\partial x_1}\right)\tag{31.59}\label{eq:trivariate_helmholtz_eq_4},\\ v_3 &= \frac{\partial\psi}{\partial x_3} + \frac{1}{\sqrt{3}}\left(\frac{\partial\chi}{\partial x_1} - \frac{\partial\chi}{\partial x_2}\right)\tag{31.60}\label{eq:trivariate_helmholtz_eq_5}. \end{align} \]
Because of
\[ \begin{align} D = \Delta\chi, & {} & \zeta = \Delta\psi \end{align} \]
the following further hold
\[ \begin{align} \Delta u_1 &= -\frac{1}{\sqrt{3}}\left(\frac{\partial\zeta}{\partial x_2} - \frac{\partial\zeta}{\partial x_3}\right) + \frac{\partial D}{\partial x_1},\tag{31.62}\label{eq:hex_laplace_u_1}\\ \Delta u_2 &= -\frac{1}{\sqrt{3}}\left(\frac{\partial\zeta}{\partial x_3} - \frac{\partial\zeta}{\partial x_1}\right) + \frac{\partial D}{\partial x_2},\tag{31.63}\label{eq:hex_laplace_u_2}\\ \Delta u_3 &= -\frac{1}{\sqrt{3}}\left(\frac{\partial\zeta}{\partial x_1} - \frac{\partial\zeta}{\partial x_2}\right) + \frac{\partial D}{\partial x_3}\tag{31.64}\label{eq:hex_laplace_u_3}, \end{align} \]
\[ \begin{align} \Delta v_1 &= \frac{\partial\zeta}{\partial x_1} + \frac{1}{\sqrt{3}}\left(\frac{\partial D}{\partial x_2} - \frac{\partial D}{\partial x_3}\right),\\ \Delta v_2 &= \frac{\partial\zeta}{\partial x_2} + \frac{1}{\sqrt{3}}\left(\frac{\partial D}{\partial x_3} - \frac{\partial D}{\partial x_1}\right),\\ \Delta v_3 &= \frac{\partial\zeta}{\partial x_3} + \frac{1}{\sqrt{3}}\left(\frac{\partial D}{\partial x_1} - \frac{\partial D}{\partial x_2}\right). \end{align} \]
Let $\alpha$ be a field of the form
\[ \begin{align} \alpha\left(\mathbf{r}\right) = \exp\left(i\mathbf{k}\cdot\mathbf{r}\right). \end{align} \]
The distance between the centers of two hexagons is denoted by $d$ (from now on this quantity is referred to as the lattice constant). One now defines, for $1 \leq j \leq 3$, the central difference quotient in the $j$-direction by
\[ \begin{align} \delta_j\alpha &= \frac{\alpha\left(\mathbf{r} + \mathbf{i}_j\frac{d}{2}\right) - \alpha\left(\mathbf{r} - \mathbf{i}_j\frac{d}{2}\right)}{d}. \end{align} \]
For this one has
\[ \begin{align} \delta_j\alpha &= \frac{\exp\left[i\mathbf{k}\cdot\left(\mathbf{r} + \mathbf{i}_j\frac{d}{2}\right)\right] - \exp\left[i\mathbf{k}\cdot\left(\mathbf{r} - \mathbf{i}_j\frac{d}{2}\right)\right]}{d}\nonumber\\ &= \exp\left(i\mathbf{k}\cdot\mathbf{r}\right)\frac{\exp\left(i\mathbf{k}\cdot\mathbf{i}_j\frac{d}{2}\right) - \exp\left(-i\mathbf{k}\cdot\mathbf{i}_j\frac{d}{2}\right)}{d} = \alpha\left(\mathbf{r}\right)\frac{\exp\left(i\mathbf{k}\cdot\mathbf{i}_j\frac{d}{2}\right) - \exp\left(-i\mathbf{k}\cdot\mathbf{i}_j\frac{d}{2}\right)}{d}\nonumber\\ &= \alpha\left(\mathbf{r}\right)\frac{\exp\left(ik_j\frac{d}{2}\right) - \exp\left(-ik_j\frac{d}{2}\right)}{d} = \alpha\left(\mathbf{r}\right)\frac{2i}{d}\sin\left(\frac{k_jd}{2}\right). \end{align} \]
If one defines
\[ \begin{align} s_j \coloneqq \sin\left(\frac{k_jd}{2}\right), \end{align} \]
one can write this in the form
\[ \begin{align} \delta_j\alpha = \frac{2i}{d}s_j\alpha \end{align} \]
One now further defines a simple averaging operator in the $j$-direction by
\[ \begin{align} \newoverline{\alpha}^{(j)} & \coloneqq \frac{\alpha\left(\mathbf{r} + \mathbf{i}_j\frac{d}{2}\right) + \alpha\left(\mathbf{r} - \mathbf{i}_j\frac{d}{2}\right)}{2} = \frac{\exp\left[i\mathbf{k}\cdot\left(\mathbf{r} + \mathbf{i}_j\frac{d}{2}\right)\right] + \exp\left[i\mathbf{k}\cdot\left(\mathbf{r} - \mathbf{i}_j\frac{d}{2}\right)\right]}{2}\nonumber\\ &= \exp\left(i\mathbf{k}\cdot\mathbf{r}\right)\frac{\exp\left(i\mathbf{k}\cdot\mathbf{i}_j\frac{d}{2}\right) + \exp\left(-i\mathbf{k}\cdot\mathbf{i}_j\frac{d}{2}\right)}{2}\nonumber\\ &= \alpha\left(\mathbf{r}\right)\frac{\exp\left(ik_j\frac{d}{2}\right) + \exp\left(-ik_j\frac{d}{2}\right)}{2} = \alpha\left(\mathbf{r}\right)\cos\left(\frac{k_jd}{2}\right).\tag{31.73}\label{eq:hex_simple_average} \end{align} \]
If one defines
\[ \begin{align} c_j \coloneqq \cos\left(\frac{k_jd}{2}\right), \end{align} \]
one can write this in the form
\[ \begin{align} \newoverline{\alpha}^{(j)} = c_j\alpha \end{align} \]
The linearized shallow water equations (13.173) - (13.174) on the f-plane can be written with the geopotential $\Phi \coloneqq gh$ in the form
\[ \begin{align} \frac{\partial\Phi}{\partial t} + \Phi_0\nabla\cdot\mathbf{v} = 0, & {} & \frac{\partial\mathbf{v}}{\partial t} + f_0\mathbf{k}\times\mathbf{v} - \nabla\Phi = 0 \end{align} \]
Here $\Phi_0 \coloneqq gD$. With Eqs. (31.33) and (31.45), this leads to the evolution equations for the components with respect to the three-element generating system:
\[ \begin{align} \frac{\partial\Phi}{\partial t} + \frac{2}{3}\Phi_0\left(\frac{\partial u_1}{\partial x_1} + \frac{\partial u_2}{\partial x_2} + \frac{\partial u_3}{\partial x_3}\right) &= 0,\\ \frac{\partial u_1}{\partial t} - \frac{f_0}{\sqrt{3}}\left(u_2 - u_3\right) + \frac{\partial\Phi}{\partial x_1} &= 0,\\ \frac{\partial u_2}{\partial t} - \frac{f_0}{\sqrt{3}}\left(u_3 - u_1\right) + \frac{\partial\Phi}{\partial x_2} &= 0,\\ \frac{\partial u_3}{\partial t} - \frac{f_0}{\sqrt{3}}\left(u_1 - u_2\right) + \frac{\partial\Phi}{\partial x_3} &= 0. \end{align} \]
Spatial discretization yields
\[ \begin{align} \frac{\partial\Phi}{\partial t} + \frac{2}{3}\Phi_0\left(\delta_1u_1 + \delta_2u_2 + \delta_3u_3\right) &= 0,\tag{31.81}\label{eq:swe_lin_hex_0}\\ \frac{\partial u_1}{\partial t} - \frac{f_0}{\sqrt{3}}\left(\newoverline{u_2}^{3} - \newoverline{u_3}^{(2)}\right) + \delta_1\Phi &= 0,\tag{31.82}\label{eq:swe_lin_hex_1}\\ \frac{\partial u_2}{\partial t} - \frac{f_0}{\sqrt{3}}\left(\newoverline{u_3}^{(1)} - \newoverline{u_1}^{3}\right) + \delta_2\Phi &= 0,\tag{31.83}\label{eq:swe_lin_hex_2}\\ \frac{\partial u_3}{\partial t} - \frac{f_0}{\sqrt{3}}\left(\newoverline{u_1}^{(2)} - \newoverline{u_2}^{(1)}\right) + \delta_3\Phi &= 0.\tag{31.84}\label{eq:swe_lin_hex_3} \end{align} \]
The arithmetic mean of the two nearest velocity components in the desired direction was used to reconstruct the wind components for computing the Coriolis acceleration. If one assumes a monochromatic plane wave for all fields, i.e.,
\[ \begin{align} \psi\left(\mathbf{r}, t\right) = \newhat{\psi}\exp\left(i\mathbf{k}\cdot\mathbf{r} - i\omega t\right) \end{align} \]
for any field $\psi$ with a complex amplitude $\newhat{\psi}$, one obtains
\[ \begin{align} -i\omega\newhat{\Phi} + \frac{4i}{3d}\Phi_0\left(s_1\newhat{u}_1 + s_2\newhat{u}_2 + s_3\newhat{u}_3\right) &= 0,\\ -i\omega\newhat{u}_1 - \frac{f_0}{\sqrt{3}}\left(c_3\newhat{u}_2 - c_2\newhat{u}_3\right) + \frac{2is_1}{d}\newhat{\Phi} &= 0,\\ -i\omega\newhat{u}_2 - \frac{f_0}{\sqrt{3}}\left(c_1\newhat{u}_3 - c_3\newhat{u}_1\right) + \frac{2is_1}{d}\newhat{\Phi} &= 0,\\ -i\omega\newhat{u}_3 - \frac{f_0}{\sqrt{3}}\left(c_2\newhat{u}_1 - c_1\newhat{u}_2\right) + \frac{2is_1}{d}\newhat{\Phi} &= 0. \end{align} \]
It should be noted that there is a phase offset between the equations, arising from the fact that the individual equations hold at different locations. This offset has been canceled out. In matrix form one obtains
\[ \begin{align} \left(\begin{array}{cccc} -i\omega & \frac{4i}{3d}\Phi_0s_1 & \frac{4i}{3d}\Phi_0s_2 & \frac{4i}{3d}\Phi_0s_3 \\ \frac{2is_1}{d} & -i\omega & -\frac{f_0}{\sqrt{3}}c_3 & \frac{f_0}{\sqrt{3}}c_2\\ \frac{2is_2}{d} & \frac{f_0}{\sqrt{3}}c_3 & -i\omega & -\frac{f_0}{\sqrt{3}}c_1\\ \frac{2is_3}{d} & -\frac{f_0}{\sqrt{3}}c_2 & \frac{f_0}{\sqrt{3}}c_1 & -i\omega \end{array}\right)\left(\begin{array}{c} \newhat{\Phi}\\ \newhat{u}_1\\ \newhat{u}_2\\ \newhat{u}_3 \end{array}\right) = \mathbf{0}. \end{align} \]
Multiplying this by the imaginary unit $i$ gives
\[ \begin{align} \left(\begin{array}{cccc} \omega & -\frac{4}{3d}\Phi_0s_1 & -\frac{4}{3d}\Phi_0s_2 & -\frac{4}{3d}\Phi_0s_3 \\ -\frac{2s_1}{d} & \omega & -\frac{if_0}{\sqrt{3}}c_3 & \frac{if_0}{\sqrt{3}}c_2\\ -\frac{2s_2}{d} & \frac{if_0}{\sqrt{3}}c_3 & \omega & -\frac{if_0}{\sqrt{3}}c_1\\ -\frac{2s_3}{d} & -\frac{if_0}{\sqrt{3}}c_2 & \frac{if_0}{\sqrt{3}}c_1 & \omega \end{array}\right)\left(\begin{array}{c} \newhat{\Phi}\\ \newhat{u}_1\\ \newhat{u}_2\\ \newhat{u}_3 \end{array}\right) = \mathbf{0}. \end{align} \]
Nontrivial solutions exist if and only if the determinant of the coefficient matrix of this system of linear equations vanishes:
\[ \begin{align} \left|\begin{array}{cccc} \omega & -\frac{4}{3d}\Phi_0s_1 & -\frac{4}{3d}\Phi_0s_2 & -\frac{4}{3d}\Phi_0s_3 \\ -\frac{2s_1}{d} & \omega & -\frac{if_0}{\sqrt{3}}c_3 & \frac{if_0}{\sqrt{3}}c_2\\ -\frac{2s_2}{d} & \frac{if_0}{\sqrt{3}}c_3 & \omega & -\frac{if_0}{\sqrt{3}}c_1\\ -\frac{2s_3}{d} & -\frac{if_0}{\sqrt{3}}c_2 & \frac{if_0}{\sqrt{3}}c_1 & \omega \end{array}\right| &\stackrel{!}{=} 0\nonumber\\ \Leftrightarrow\omega\left|\begin{array}{ccc} \omega & -\frac{if_0}{\sqrt{3}}c_3 & \frac{if_0}{\sqrt{3}}c_2\\ \frac{if_0}{\sqrt{3}}c_3 & \omega & -\frac{if_0}{\sqrt{3}}c_1\\ -\frac{if_0}{\sqrt{3}}c_2 & \frac{if_0}{\sqrt{3}}c_1 & \omega \end{array} \right| + \frac{4}{3d}\Phi_0s_1\left|\begin{array}{ccc} -\frac{2s_1}{d} & -\frac{if_0}{\sqrt{3}}c_3 & \frac{if_0}{\sqrt{3}}c_2\\ -\frac{2s_2}{d} & \omega & -\frac{if_0}{\sqrt{3}}c_1\\ -\frac{2s_3}{d} & \frac{if_0}{\sqrt{3}}c_1 & \omega \end{array}\right| & \nonumber\\ - \frac{4}{3d}\Phi_0s_2\left|\begin{array}{ccc} -\frac{2s_1}{d} & \omega & \frac{if_0}{\sqrt{3}}c_2\\ -\frac{2s_2}{d} & \frac{if_0}{\sqrt{3}}c_3 & -\frac{if_0}{\sqrt{3}}c_1\\ -\frac{2s_3}{d} & -\frac{if_0}{\sqrt{3}}c_2 & \omega \end{array} \right| + \frac{4}{3d}\Phi_0s_3\left|\begin{array}{ccc} -\frac{2s_1}{d} & \omega & -\frac{if_0}{\sqrt{3}}c_3\\ -\frac{2s_2}{d} & \frac{if_0}{\sqrt{3}}c_3 & \omega\\ -\frac{2s_3}{d} & -\frac{if_0}{\sqrt{3}}c_2 & \frac{if_0}{\sqrt{3}}c_1 \end{array} \right| &= 0 \end{align} \] \[ \begin{align} \Leftrightarrow \omega^4 - \frac{\omega^2f_0^2}{3}\left(c_1^2 + c_2^2 + c_3^2\right) + \frac{4}{3d}\Phi_0s_1\left|\begin{array}{ccc} -\frac{2s_1}{d} & -\frac{if_0}{\sqrt{3}}c_3 & \frac{if_0}{\sqrt{3}}c_2\\ -\frac{2s_2}{d} & \omega & -\frac{if_0}{\sqrt{3}}c_1\\ -\frac{2s_3}{d} & \frac{if_0}{\sqrt{3}}c_1 & \omega \end{array}\right| & \nonumber\\ -\frac{4}{3d}\Phi_0s_2\left|\begin{array}{ccc} -\frac{2s_1}{d} & \omega & \frac{if_0}{\sqrt{3}}c_2\\ -\frac{2s_2}{d} & \frac{if_0}{\sqrt{3}}c_3 & -\frac{if_0}{\sqrt{3}}c_1\\ -\frac{2s_3}{d} & -\frac{if_0}{\sqrt{3}}c_2 & \omega \end{array}\right| + \frac{4}{3d}\Phi_0s_3\left|\begin{array}{ccc} -\frac{2s_1}{d} & \omega & -\frac{if_0}{\sqrt{3}}c_3\\ -\frac{2s_2}{d} & \frac{if_0}{\sqrt{3}}c_3 & \omega\\ -\frac{2s_3}{d} & -\frac{if_0}{\sqrt{3}}c_2 & \frac{if_0}{\sqrt{3}}c_1 \end{array} \right| &= 0 \end{align} \]
One has
\[ \begin{align} \frac{4}{3d}\Phi_0s_1\left|\begin{array}{ccc} -\frac{2s_1}{d} & -\frac{if_0}{\sqrt{3}}c_3 & \frac{if_0}{\sqrt{3}}c_2\\ -\frac{2s_2}{d} & \omega & -\frac{if_0}{\sqrt{3}}c_1\\ -\frac{2s_3}{d} & \frac{if_0}{\sqrt{3}}c_1 & \omega \end{array}\right| &= -\frac{8\omega^2}{3d^2}\Phi_0s_1^2 + \frac{8\Phi_0f_0^2}{9d^2}\left(s_1c_1s_3c_3 + s_1c_1s_2c_2 + s_1c_1s_1c_1\right) + \frac{8if_0}{3\sqrt{3}d}s_1\left(s_3c_2 - c_3s_2\right)\nonumber\\ &= -\frac{8\omega^2}{3d^2}\Phi_0s_1^2 + \frac{8\Phi_0f_0^2}{9d^2}s_1c_1\left(s_1c_1 + s_2c_2 + s_3c_3\right) - \frac{8if_0}{3\sqrt{3}d}s_1\left(s_2c_3 - s_3c_2\right). \end{align} \]
The remaining two partial determinants arise analogously by cyclically exchanging the indices:
\[ \begin{align} -\frac{4}{3d}\Phi_0s_2\left|\begin{array}{ccc} -\frac{2s_1}{d} & \omega & \frac{if_0}{\sqrt{3}}c_2\\ -\frac{2s_2}{d} & \frac{if_0}{\sqrt{3}}c_3 & -\frac{if_0}{\sqrt{3}}c_1\\ -\frac{2s_3}{d} & -\frac{if_0}{\sqrt{3}}c_2 & \omega \end{array}\right| &= -\frac{8\omega^2}{3d^2}\Phi_0s_2^2 + \frac{8\Phi_0f_0^2}{9d^2}s_2c_2\left(s_1c_1 + s_2c_2 + s_3c_3\right) - \frac{8if_0}{3\sqrt{3}d}s_2\left(s_3c_1 - s_1c_3\right),\nonumber\\ \frac{4}{3d}\Phi_0s_3\left|\begin{array}{ccc} -\frac{2s_1}{d} & \omega & -\frac{if_0}{\sqrt{3}}c_3\\ -\frac{2s_2}{d} & \frac{if_0}{\sqrt{3}}c_3 & \omega\\ -\frac{2s_3}{d} & -\frac{if_0}{\sqrt{3}}c_2 & \frac{if_0}{\sqrt{3}}c_1 \end{array} \right| &= -\frac{8\omega^2}{3d^2}\Phi_0s_3^2 + \frac{8\Phi_0f_0^2}{9d^2}s_3c_3\left(s_1c_1 + s_2c_2 + s_3c_3\right) - \frac{8if_0}{3\sqrt{3}d}s_3\left(s_1c_2 - s_2c_1\right). \end{align} \]
Thus, the characteristic polynomial finally reads:
For the case of well-resolved waves $\left|\mathbf{k}\right| \ll \frac{1}{d}$, one has
\[ \begin{align} c_j \approx 1, & {} & s_j \approx \frac{k_jd}{2}. \end{align} \]
Substituting this into Eq. (31.96), one obtains
\[ \begin{align} \omega^4 - \omega^2\left[f_0^2 + \frac{8\Phi_0}{3d^2}\left(\frac{k_1^2d^2}{4} + \frac{k_2^2d^2}{4} + \frac{k_3^2d^2}{4}\right)\right] + \frac{8\Phi_0f_0^2}{9d^2}\underbrace{\left(\frac{k_1d}{2} + \frac{k_1d}{2} + \frac{k_1d}{2}\right)^2}_{\approx 0} &\approx 0\nonumber\\ \Leftrightarrow\omega^4 - \omega^2\left[f_0^2 + \frac{2\Phi_0}{3}\left(k_1^2 + k_2^2 + k_3^2\right)\right] &\approx 0.\nonumber\\ \Leftrightarrow\omega^4 - \omega^2\left[f_0^2 + \Phi\mathbf{k}^2\right] &\approx 0. \end{align} \]
The last step took advantage of the fact that $k_1^2 + k_2^2 + k_3^2 = \frac{3}{2}\mathbf{k}^2$. This can be seen by assuming, without loss of generality, that $\mathbf{k} \parallel \mathbf{i}_1$, since then $k_1 = k$, $k_2 = -k\sin\left(30^\circ\right) = -\frac{k}{2}$ and $k_3 = -\frac{k}{2}$. This leads to the two inertia-gravity modes
\[ \begin{align} \omega^2 \approx f_0^2 + \Phi_0\left|\mathbf{k}\right|^2 \end{align} \]
as well as the geostrophic mode
\[ \begin{align} \omega \approx 0. \end{align} \]
However, for less well-resolved waves, the term $\propto\left(s_1c_1 + s_2c_2 + s_3c_3\right)^2$ causes the geostrophic mode to split into two submodes with $\omega \not= 0$. This is problematic and corresponds to the findings of Sect. 30.9.
In a planar triangular or hexagonal C-grid discretization, the velocity components are rotated by 120 or 60$^\circ$ with respect to each other, but the components are not defined at the same location. Therefore, the condition Eq. (31.18) does not hold straightforwardly for the vector components, but only after an interpolation to reference locations, for which the centers of the cells are used here. This interpolation is done using operators
\[ \begin{align} \newtilde{u}^{(i)} \end{align} \]
with $1 \leq i \leq 3$, where $u$ is any vector component. A separate averaging operator is permitted here for each spatial direction $i$. The conditions Eq. (31.18) can be expressed here in the form
\[ \begin{align} \newtilde{u}_1^{(1)} + \newtilde{u}_2^{(2)} + \newtilde{u}_3^{(3)} \hastobe 0, & {} & \newtilde{v}_1^{(1)} + \newtilde{v}_2^{(2)} + \newtilde{v}_3^{3} \hastobe 0\tag{31.102}\label{eq:linear_condition_trivariate_discrete} \end{align} \]
These equations are to be understood as requirements for the vector fields $\mathbf{v}$. If gradient fields $\nabla\alpha$ do not satisfy the condition Eq. (31.102), one has a problem, since the computation of gradient fields on the C-grid is fixed. Therefore, the averaging of Eq. (31.22) is used as a requirement for the averaging operator to be found:
For the Helmholtz decomposition, Eqs. (31.55) - (31.60), averaging in the respectively relevant spatial directions yields
\[ \begin{align} \newtilde{u}_1^{(1)} &= \frac{1}{\sqrt{3}}\left(\newtilde{\delta_3\psi}^{(1)} - \newtilde{\delta_2\psi}^{(1)}\right) + \newtilde{\delta_1\chi}^{(1)},\\ \newtilde{u}_2^{(2)} &= \frac{1}{\sqrt{3}}\left(\newtilde{\delta_1\psi}^{(2)} - \newtilde{\delta_3\psi}^{(2)}\right) + \newtilde{\delta_2\chi}^{(2)},\\ \newtilde{u}_3^{3} &= \frac{1}{\sqrt{3}}\left(\newtilde{\delta_2\psi}^{3} - \newtilde{\delta_1\psi}^{3}\right) + \newtilde{\delta_3\chi}^{3},\\ \newtilde{v}_1^{(1)} &= \frac{1}{\sqrt{3}}\left(\newtilde{\delta_2\chi}^{(1)} - \newtilde{\delta_3\chi}^{(1)}\right) + \newtilde{\delta_1\psi}^{(1)},\\ \newtilde{v}_2^{(2)} &= \frac{1}{\sqrt{3}}\left(\newtilde{\delta_3\chi}^{(2)} - \newtilde{\delta_1\chi}^{(2)}\right) + \newtilde{\delta_2\psi}^{(2)},\\ \newtilde{v}_3^{3} &= \frac{1}{\sqrt{3}}\left(\newtilde{\delta_1\chi}^{3} - \newtilde{\delta_2\chi}^{3}\right) + \newtilde{\delta_3\psi}^{3}. \end{align} \]
This does not yet satisfy Eqs. (31.102) - (31.102); rather, a preparatory averaging of the scalar fields is required:
\[ \begin{align} \newtilde{u}_1^{(1)} &= \frac{1}{\sqrt{3}}\left(\newtilde{\delta_3\newtilde{\psi}^{(2)}}^{(1)} - \newtilde{\delta_2\newtilde{\psi}^{3}}^{(1)}\right) + \newtilde{\delta_1\chi}^{(1)}\\ \newtilde{u}_2^{(2)} &= \frac{1}{\sqrt{3}}\left(\newtilde{\delta_1\newtilde{\psi}^{3}}^{(2)} - \newtilde{\delta_3\newtilde{\psi}^{(1)}}^{(2)}\right) + \newtilde{\delta_2\chi}^{(2)}\\ \newtilde{u}_3^{3} &= \frac{1}{\sqrt{3}}\left(\newtilde{\delta_2\newtilde{\psi}^{(1)}}^{3} - \newtilde{\delta_1\newtilde{\psi}^{(2)}}^{3}\right) + \newtilde{\delta_3\chi}^{3}\\ \newtilde{v}_1^{(1)} &= \frac{1}{\sqrt{3}}\left(\newtilde{\delta_2\newtilde{\chi}^{3}}^{(1)} - \newtilde{\delta_3\newtilde{\chi}^{(2)}}^{(1)}\right) + \newtilde{\delta_1\psi}^{(1)}\\ \newtilde{v}_2^{(2)} &= \frac{1}{\sqrt{3}}\left(\newtilde{\delta_3\newtilde{\chi}^{(1)}}^{(2)} - \newtilde{\delta_1\newtilde{\chi}^{3}}^{(2)}\right) + \newtilde{\delta_2\psi}^{(2)}\\ \newtilde{v}_3^{3} &= \frac{1}{\sqrt{3}}\left(\newtilde{\delta_1\newtilde{\chi}^{(2)}}^{3} - \newtilde{\delta_2\newtilde{\chi}^{(1)}}^{3}\right) + \newtilde{\delta_3\psi}^{3} \end{align} \]
This implies
\[ \begin{align} u_1 &= \frac{1}{\sqrt{3}}\left(\delta_3\newtilde{\psi}^{(2)} - \delta_2\newtilde{\psi}^{3}\right) + \delta_1\chi\tag{31.116}\label{eq:hex:helmholtz_0},\\ u_2 &= \frac{1}{\sqrt{3}}\left(\delta_1\newtilde{\psi}^{3} - \delta_3\newtilde{\psi}^{(1)}\right) + \delta_2\chi\tag{31.117}\label{eq:hex:helmholtz_1},\\ u_3 &= \frac{1}{\sqrt{3}}\left(\delta_2\newtilde{\psi}^{(1)} - \delta_1\newtilde{\psi}^{(2)}\right) + \delta_3\chi\tag{31.118}\label{eq:hex:helmholtz_2},\\ v_1 &= \frac{1}{\sqrt{3}}\left(\delta_2\newtilde{\chi}^{3} - \delta_3\newtilde{\chi}^{(2)}\right) + \delta_1\psi,\\ v_2 &= \frac{1}{\sqrt{3}}\left(\delta_3\newtilde{\chi}^{(1)} - \delta_1\newtilde{\chi}^{3}\right) + \delta_2\psi,\\ v_3 &= \frac{1}{\sqrt{3}}\left(\delta_1\newtilde{\chi}^{(2)} - \delta_2\newtilde{\chi}^{(1)}\right) + \delta_3\psi. \end{align} \]
Now the question arises what the averaging operator $\newtilde{\psi}^{(j)}$ must concretely look like. For the operator introduced in Eq. (31.73), one has
\[ \begin{align} &\newoverline{\delta_1\psi}^{(1)} + \newoverline{\delta_1\psi}^{(2)} + \newoverline{\delta_1\psi}^{3}\nonumber\\ &= \psi\left(\mathbf{r}\right)\frac{2i}{d}\left[\sin\left(\frac{k_1d}{2}\right)\cos\left(\frac{k_1d}{2}\right) + \sin\left(\frac{k_2d}{2}\right)\cos\left(\frac{k_2d}{2}\right) + \sin\left(\frac{k_3d}{2}\right)\cos\left(\frac{k_3d}{2}\right)\right]\nonumber\\ &= \psi\left(\mathbf{r}\right)\frac{i}{d}\left[2\sin\left(\frac{k_1d}{2}\right)\cos\left(\frac{k_1d}{2}\right) + 2\sin\left(\frac{k_2d}{2}\right)\cos\left(\frac{k_2d}{2}\right) + 2\sin\left(\frac{k_3d}{2}\right)\cos\left(\frac{k_3d}{2}\right)\right]\nonumber\\ &= \psi\left(\mathbf{r}\right)\frac{i}{d}\left[\sin\left(k_1d\right) + \sin\left(k_2d\right) + \sin\left(k_3d\right)\right] \stackrel{\text{in general}}{\not=} 0. \end{align} \]
This operator therefore does not satisfy the condition Eq. (31.103). Therefore, a second averaging operator is defined
\[ \begin{align} \newoverline{\psi}^{((1))} & \coloneqq \newoverline{\psi}^{(2, 3)} = \newoverline{\newoverline{\psi}^{3}}^{(2)} = \newoverline{\newoverline{\psi}^{(2)}}^{3} = \psi\left(\mathbf{r}\right)\cos\left(\frac{k_2d}{2}\right)\cos\left(\frac{k_3d}{2}\right).\text{ (and cyclic)} \end{align} \]
For this one has
\[ \begin{align} \newoverline{\delta_1\alpha}^{((1))} + \newoverline{\delta_2\alpha}^{((2))} + \newoverline{\delta_3\alpha}^{(3)} &= \alpha\left(\mathbf{r}\right)\frac{2i}{d}\Big[\sin\left(\frac{k_1d}{2}\right)\cos\left(\frac{k_2d}{2}\right)\cos\left(\frac{k_3d}{2}\right) + \cos\left(\frac{k_1d}{2}\right)\sin\left(\frac{k_2d}{2}\right)\cos\left(\frac{k_3d}{2}\right)\nonumber\\ & + \cos\left(\frac{k_1d}{2}\right)\cos\left(\frac{k_2d}{2}\right)\sin\left(\frac{k_3d}{2}\right)\Big]\nonumber\\ &= \alpha\left(\mathbf{r}\right)\frac{2i}{d}\Big[\sin\left(\frac{k_1d + k_2d}{2}\right)\cos\left(\frac{k_3d}{2}\right) + \cos\left(\frac{k_1d}{2}\right)\cos\left(\frac{k_2d}{2}\right)\sin\left(\frac{k_3d}{2}\right)\Big]\nonumber\\ &= \alpha\left(\mathbf{r}\right)\frac{2i}{d}\Big[\sin\left(\frac{k_1 + k_2 + k_3}{2}d\right) - \cos\left(\frac{k_1 + k_2}{2}d\right)\sin\left(\frac{k_3}{2}d\right)\nonumber\\ &+ \cos\left(\frac{k_1d}{2}\right)\cos\left(\frac{k_2d}{2}\right)\sin\left(\frac{k_3d}{2}\right)\Big]. \end{align} \]
With $k_1 + k_2 + k_3 = 0$ follows
\[ \begin{align} \newoverline{\delta_1\alpha}^{((1))} + \newoverline{\delta_2\alpha}^{((2))} + \newoverline{\delta_3\alpha}^{(3)} &= \alpha\left(\mathbf{r}\right)\frac{2i}{d}\Big[-\cos\left(\frac{k_1 + k_2}{2}d\right)\sin\left(\frac{k_3}{2}d\right) + \cos\left(\frac{k_1d}{2}\right)\cos\left(\frac{k_2d}{2}\right)\sin\left(\frac{k_3d}{2}\right)\Big]\nonumber\\ &= \alpha\left(\mathbf{r}\right)\frac{2i}{d}\sin\left(\frac{k_3}{2}d\right)\Big[\cos\left(\frac{k_1d}{2}\right)\cos\left(\frac{k_2d}{2}\right) - \cos\left(\frac{k_1 + k_2}{2}d\right)\Big]\nonumber\\ &= \alpha\left(\mathbf{r}\right)\frac{2i}{d}\sin\left(\frac{k_3}{2}d\right)\Big[\cos\left(\frac{k_1d}{2}\right)\cos\left(\frac{k_2d}{2}\right) - \cos\left(\frac{k_1d}{2}\right)\cos\left(\frac{k_2d}{2}\right)\nonumber\\ &+ \sin\left(\frac{k_1d}{2}\right)\sin\left(\frac{k_2d}{2}\right)\Big]\nonumber\\ &= \alpha\left(\mathbf{r}\right)\frac{2i}{d}\sin\left(\frac{k_1}{2}d\right)\sin\left(\frac{k_2d}{2}\right)\sin\left(\frac{k_3d}{2}\right) \stackrel{\text{in general}}{\not=} 0.\tag{31.125}\label{eq:thuburn_op_deriv_0} \end{align} \]
This operator also does not satisfy the condition Eq. (31.103). However, one continues the calculation, again using $k_1 + k_2 + k_3 = 0$:
\[ \begin{align} & \newoverline{\delta_1\alpha}^{((1))} + \newoverline{\delta_2\alpha}^{((2))} + \newoverline{\delta_3\alpha}^{(3)} = \alpha\left(\mathbf{r}\right)\frac{2i}{d}\sin\left(\frac{k_1}{2}d\right)\sin\left(\frac{k_2d}{2}\right)\sin\left(\frac{k_3d}{2}\right)\nonumber\\ &= \alpha\left(\mathbf{r}\right)\frac{2i}{d}\sin\left(\frac{k_1}{2}d\right)\sin\left(\frac{k_2d}{2}\right)\sin\left(-\frac{k_1 + k_2}{2}d\right) = -\alpha\left(\mathbf{r}\right)\frac{2i}{d}\sin\left(\frac{k_1}{2}d\right)\sin\left(\frac{k_2d}{2}\right)\sin\left(\frac{k_1 + k_2}{2}d\right)\nonumber\\ &= -\alpha\left(\mathbf{r}\right)\frac{2i}{d}\sin\left(\frac{k_1}{2}d\right)\sin\left(\frac{k_2d}{2}\right)\left[\sin\left(\frac{k_1}{2}d\right)\cos\left(\frac{k_2}{2}d\right) + \cos\left(\frac{k_1}{2}d\right)\sin\left(\frac{k_2}{2}d\right)\right]\nonumber\\ &= -\alpha\left(\mathbf{r}\right)\frac{2i}{d}\left[\sin\left(\frac{k_2}{2}d\right)\sin\left(\frac{k_1}{2}d\right)^2\cos\left(\frac{k_2}{2}d\right) + \sin\left(\frac{k_1}{2}d\right)\cos\left(\frac{k_1}{2}d\right)\sin\left(\frac{k_2}{2}d\right)^2\right]\nonumber\\ &= -\alpha\left(\mathbf{r}\right)\frac{i}{d}\left[\sin\left(k_2d\right)\sin\left(\frac{k_1}{2}d\right)^2 + \sin\left(k_1d\right)\sin\left(\frac{k_2}{2}d\right)^2\right]\nonumber\\ &= -\alpha\left(\mathbf{r}\right)\frac{i}{d}\left[\sin\left(k_2d\right) + \sin\left(k_1d\right) - \sin\left(k_2d\right)\cos\left(\frac{k_1}{2}d\right)^2 - \sin\left(k_1d\right)\cos\left(\frac{k_2}{2}d\right)^2\right]\nonumber \end{align} \] \[ \begin{align} &= -\alpha\left(\mathbf{r}\right)\frac{i}{d}\Big[\sin\left(k_2d\right) + \sin\left(k_1d\right) - \sin\left(k_2d\right)\cos\left(k_1\right) - \sin\left(k_1d\right)\cos\left(k_2d\right) + \sin\left(k_2d\right)\sin\left(\frac{k_1}{2}d\right)^2\nonumber\\ & + \sin\left(k_1d\right)\sin\left(\frac{k_2}{2}d\right)^2\Big]\nonumber\\ &= -\alpha\left(\mathbf{r}\right)\frac{i}{d}\Big[\sin\left(k_2d\right) + \sin\left(k_1d\right) - \sin\left[\left(k_1 + k_2\right)d\right] + \sin\left(k_2d\right)\sin\left(\frac{k_1}{2}d\right)^2 + \sin\left(k_1d\right)\sin\left(\frac{k_2}{2}d\right)^2\Big]\nonumber\\ &= -\alpha\left(\mathbf{r}\right)\frac{i}{d}\Big[\sin\left(k_2d\right) + \sin\left(k_1d\right) - \sin\left(-k_3d\right) + \sin\left(k_2d\right)\sin\left(\frac{k_1}{2}d\right)^2 + \sin\left(k_1d\right)\sin\left(\frac{k_2}{2}d\right)^2\Big]\nonumber\\ &= -\alpha\left(\mathbf{r}\right)\frac{i}{d}\Big[\sin\left(k_2d\right) + \sin\left(k_1d\right) + \sin\left(k_3d\right) + \sin\left(k_2d\right)\sin\left(\frac{k_1}{2}d\right)^2 + \sin\left(k_1d\right)\sin\left(\frac{k_2}{2}d\right)^2\Big]\nonumber\\ &= -\left(\newoverline{\delta_1\alpha}^{(1)} + \newoverline{\delta_1\alpha}^{(2)} + \newoverline{\delta_1\alpha}^{3}\right) - \alpha\left(\mathbf{r}\right)\frac{i}{d}\Big[\sin\left(k_2d\right)\sin\left(\frac{k_1}{2}d\right)^2 + \sin\left(k_1d\right)\sin\left(\frac{k_2}{2}d\right)^2\Big] \end{align} \]
One now computes
\[ \begin{align} &\sin\left(k_2d\right)\sin\left(\frac{k_1}{2}d\right)^2 + \sin\left(k_1d\right)\sin\left(\frac{k_2}{2}d\right)^2 = \sin\left(\frac{k_2}{2}d + \frac{k_2}{2}d\right)\sin\left(\frac{k_1}{2}d\right)^2 + \sin\left(\frac{k_1}{2}d + \frac{k_1}{2}d\right)\sin\left(\frac{k_2}{2}d\right)^2\nonumber\\ &= 2\sin\left(\frac{k_1}{2}d\right)^2\sin\left(\frac{k_2}{2}d\right)\cos\left(\frac{k_2}{2}d\right) + 2\sin\left(\frac{k_2}{2}d\right)^2\sin\left(\frac{k_1}{2}d\right)\cos\left(\frac{k_1}{2}d\right)\nonumber\\ &= \sin\left(\frac{k_1}{2}d\right)\sin\left(\frac{k_2}{2}d\right)\left[2\sin\left(\frac{k_1}{2}d\right)\cos\left(\frac{k_2}{2}d\right) + 2\cos\left(\frac{k_1}{2}d\right)\sin\left(\frac{k_2}{2}d\right)\right]\nonumber\\ &= \sin\left(\frac{k_1}{2}d\right)\sin\left(\frac{k_2}{2}d\right)\sin\left(\frac{k_1 + k_2}{2}d\right) = \sin\left(\frac{k_1}{2}d\right)\sin\left(\frac{k_2}{2}d\right)\sin\left(-\frac{k_3}{2}d\right)\nonumber\\ &= -\sin\left(\frac{k_1}{2}d\right)\sin\left(\frac{k_2}{2}d\right)\sin\left(\frac{k_3}{2}d\right). \end{align} \]
Thus follows
\[ \begin{align} & \newoverline{\delta_1\alpha}^{((1))} + \newoverline{\delta_2\alpha}^{((2))} + \newoverline{\delta_3\alpha}^{(3)} = -\left(\newoverline{\delta_1\alpha}^{(1)} + \newoverline{\delta_1\alpha}^{(2)} + \newoverline{\delta_1\alpha}^{3}\right) + \alpha\left(\mathbf{r}\right)\frac{i}{d}\sin\left(\frac{k_1}{2}d\right)\sin\left(\frac{k_2}{2}d\right)\sin\left(\frac{k_3}{2}d\right)\nonumber\\ &\stackrel{\href{#eq:thuburn_op_deriv_0}{\text{Eq. (31.125)}}}{=} -\left(\newoverline{\delta_1\alpha}^{(1)} + \newoverline{\delta_1\alpha}^{(2)} + \newoverline{\delta_1\alpha}^{3}\right) + \frac{1}{2}\left(\newoverline{\delta_1\alpha}^{((1))} + \newoverline{\delta_2\alpha}^{((2))} + \newoverline{\delta_3\alpha}^{(3)}\right)\nonumber\\ &\Rightarrow\newoverline{\delta_1\alpha}^{((1))} + \newoverline{\delta_2\alpha}^{((2))} + \newoverline{\delta_3\alpha}^{(3)} = -\frac{1}{2}\left(\newoverline{\delta_1\alpha}^{(1)} + \newoverline{\delta_1\alpha}^{(2)} + \newoverline{\delta_1\alpha}^{3}\right).\tag{31.128}\label{eq:thuburn_op_deriv_1} \end{align} \]
One now makes the ansatz
\[ \begin{align} \newtilde{\alpha}^{(j)} = \beta\newoverline{\alpha}^{(j)} + \left(1 - \beta\right)\newoverline{\alpha}^{((j))} \end{align} \]
with $0 \leq \beta \leq 0$. Substituting Eq. (31.128) here, one obtains
\[ \begin{align} \newtilde{\delta_1\alpha}^{(1)} + \newtilde{\delta_2\alpha}^{(2)} + \newtilde{\delta_3\alpha}^{3} = \left(\beta - \frac{1}{2}\left(1 - \beta\right)\right)\left(\newoverline{\delta_1\alpha}^{(1)} + \newoverline{\delta_1\alpha}^{(2)} + \newoverline{\delta_1\alpha}^{3}\right) \stackrel{!}{=} 0. \end{align} \]
This is fulfilled for
\[ \begin{align} \beta - \frac{1}{2}\left(1 - \beta\right) &= 0 \Leftrightarrow \beta - \frac{1}{2} + \frac{\beta}{2} = \frac{3\beta}{2} - \frac{1}{2} = 0\nonumber\\ \Leftrightarrow \beta = \frac{1}{3}. \end{align} \]
Thus, one has
This operator is called the Thuburn operator because it was first used by Thuburn in [38].
Applying the Thuburn operator in the linearized shallow water equations on the hexagonal grid, Eqs. (31.81) - (31.84), leads to
\[ \begin{align} \frac{\partial\Phi}{\partial t} + \frac{2}{3}\Phi_0\left(\delta_1u_1 + \delta_2u_2 + \delta_3u_3\right) &= 0,\\ \frac{\partial u_1}{\partial t} - \frac{f_0}{\sqrt{3}}\left(\newtilde{u}_2^{(3)} - \newtilde{u}_3^{(2)}\right) + \delta_1\Phi &= 0,\tag{31.134}\label{eq:x_1_momentum_thuburn_c-grid}\\ \frac{\partial u_2}{\partial t} - \frac{f_0}{\sqrt{3}}\left(\newtilde{u}_3^{(1)} - \newtilde{u}_1^{(3)}\right) + \delta_2\Phi &= 0,\\ \frac{\partial u_3}{\partial t} - \frac{f_0}{\sqrt{3}}\left(\newtilde{u}_1^{(2)} - \newtilde{u}_2^{(1)}\right) + \delta_3\Phi &= 0. \end{align} \]
From this it follows immediately
\[ \begin{align} \frac{\partial}{\partial t}\left(\newtilde{u}_1^{(1)} + \newtilde{u}_2^{(2)} + \newtilde{u}_3^{(3)}\right) = 0.\tag{31.137}\label{eq:thuburn_con_swe_c-grid} \end{align} \]
The dispersion relation of this system of equations can be calculated analogously to the derivation of Eq. (31.96) by performing the substitution
\[ \begin{align} c_1 \to a_1 \coloneqq\frac{c_1 + 2c_2c_3}{3}\text{ (and cyclic)} \end{align} \]
This leads to
\[ \begin{align} \omega^4 - \omega^2\left[\frac{f_0^2}{3}\left(a_1^2 + a_2^2 + a_3^2\right) + \frac{8\Phi_0}{3d^2}\left(s_1^2 + s_2^2 + s_3^2\right)\right] + \frac{8\Phi_0f_0^2}{9d^2}\left(s_1a_1 + s_2a_2 + s_3a_3\right)^2 = 0.\tag{31.139}\label{eq:disp_rel_hex_c_mod_0} \end{align} \]
If one sets $\alpha = \exp\left(i\mathbf{k}\cdot\mathbf{r}\right)$ in Eq. (31.103), one obtains
\[ \begin{align} \frac{2i}{d}\left(a_1s_1 + a_2s_2 + a_3s_3\right)\alpha = 0 \Leftrightarrow a_1s_1 + a_2s_2 + a_3s_3 = 0. \end{align} \]
Thus, Eq. (31.139) becomes
\[ \begin{align} \omega^4 - \omega^2\left[\frac{f_0^2}{3}\left(a_1^2 + a_2^2 + a_3^2\right) + \frac{8\Phi_0}{3d^2}\left(s_1^2 + s_2^2 + s_3^2\right)\right] = 0.\tag{31.141}\label{eq:thuburn_leads_to_geostrophic_mode} \end{align} \]
If one uses Eq. (31.132) to reconstruct the tangential velocity components, a stationary geostrophic mode thus arises.
On a regular hexagonal lattice with lattice constant $d$, a triangle has area
\[ \begin{align} A_v = \frac{d}{2}l. \end{align} \]
Here one has
\[ \begin{align} l = \cos\left(30^\circ\right)d = \frac{\sqrt{3}}{2}d, \end{align} \]
from which
\[ \begin{align} A_v = \frac{d^2\sqrt{3}}{4}\tag{31.144}\label{eq:hex_face_triangle} \end{align} \]
follows. Thus, on a triangle (superscript $t$) pointing upward (index $u$), one can compute the divergence according to Gauss's theorem in the form
\[ \begin{align} D^{(t)}_u = -\frac{dv_1 + dv_2 + dv_3}{A_v} = -\frac{4}{\sqrt{3}d}\left(v_1 + v_2 + v_3\right) \end{align} \]
The minus sign arises from the fact that in such a triangle the coordinate axes of the $j$-system point inward. For a triangle pointing downward (index $l$), it is the other way around; one has
\[ \begin{align} D^{(t)}_l = \frac{dv_1 + dv_2 + dv_3}{A_v} = \frac{4}{\sqrt{3}d}\left(v_1 + v_2 + v_3\right). \end{align} \]
This can be summarized in the form
\[ \begin{align} D^{(t)}_{u, l} = \mp\frac{dv_1 + dv_2 + dv_3}{A_v} = \mp\frac{4}{\sqrt{3}d}\left(v_1 + v_2 + v_3\right) = \mp\frac{4}{\sqrt{3}d}\sum_{e \in t}v_e \end{align} \]
where $e$ denotes the edges of the triangle. Analogously, for the relative vorticity, Stokes' theorem gives
\[ \begin{align} \zeta^{(t)}_{u, l} = \pm\frac{du_1 + du_2 + du_3}{A_v} = \pm\frac{4}{\sqrt{3}d}\sum_{e \in t}u_e. \end{align} \]
Let $h$ now denote a hexagon. Then $u \in h$ denote the three upward-pointing triangles that overlap with the hexagon $h$ under consideration and similarly $l \in h$ denote the three downward-pointing triangles that overlap with $h$. With these notations one can write
\[ \begin{align} \frac{\sqrt{3}d}{4}\left(\sum_{l \in h}D_l^{(t)} - \sum_{u \in h}D_u^{(t)}\right) &= \sum_{l \in h}\sum_{e \in l}v_e + \sum_{u \in h}\sum_{e \in u}v_e = \sum_{l \in h}\left(\sum_{e \in l, e \in h}v_e + \sum_{e \in l, e \not\in h}v_e\right) + \sum_{u \in h}\left(\sum_{e \in u, e \in h}v_e + \sum_{e \in u, e \not\in h}v_e\right) \end{align} \]
In the last transformation step, the edges of the triangles were divided into those that are also edges of the hexagon under consideration ($e \in h$) and those that are not ($e \not\in h$). Factoring out yields
\[ \begin{align} \frac{\sqrt{3}d}{4}\left(\sum_{l \in h}D_l^{(t)} - \sum_{u \in h}D_u^{(t)}\right) &= \sum_{l \in h}\sum_{e \in l, e \in h}v_e + \sum_{l \in h}\sum_{e \in l, e \not\in h}v_e + \sum_{u \in h}\sum_{e \in u, e \in h}v_e + \sum_{u \in h}\sum_{e \in u, e \not\in h}v_e. \end{align} \]
As one can see, all velocity components contribute with a positive sign to $\sum_{l \in h}D_l^{(t)} - \sum_{u \in h}D_u^{(t)}$. In the sums $\sum_{l \in h}\sum_{e \in l, e \in h}$ and $\sum_{u \in h}\sum_{e \in u, e \in h}$, each edge of the hexagon $h$ occurs exactly once. Thus one has
\[ \begin{align} \frac{\sqrt{3}d}{4}\left(\sum_{l \in h}D_l^{(t)} - \sum_{u \in h}D_u^{(t)}\right) &= 2\sum_{e \in h}v_e + \sum_{l \in h}\sum_{e \in l, e \not\in h}v_e + \sum_{u \in h}\sum_{e \in u, e \not\in h}v_e. \end{align} \]
Those edges that do not border the hexagon under consideration each contribute exactly once to $\sum_{l \in h}D_l^{(t)} - \sum_{u \in h}D_u^{(t)}$:
\[ \begin{align} \frac{\sqrt{3}d}{4}\left(\sum_{l \in h}D_l^{(t)} - \sum_{u \in h}D_u^{(t)}\right) &= 2\sum_{e \in h}v_e + \sum_{t \in h}\sum_{e \in t, e \not\in h}v_e \end{align} \]
At this point, consider the quantity $\newtilde{v}_1^{(1)}$, which is defined at the center of the hexagon under consideration. It is a superposition of four velocity components, two of which lie on the $\mathbf{i}_1$-axis and two perpendicular to it. The former contribute to $\newtilde{v}_1^{(1)}$ with weight $1/3$, the latter with weight $1/6.$ This can also be transferred cyclically to the other spatial directions. Thus one has
\[ \begin{align} \frac{\sqrt{3}d}{4}\left(\sum_{l \in h}D_l^{(t)} - \sum_{u \in h}D_u^{(t)}\right) &= 6\left(\newtilde{v}_1^{(1)} + \newtilde{v}_2^{(2)} + \newtilde{v}_3^{3}\right).\tag{31.153}\label{eq:checkerboard_triangular} \end{align} \]
If the condition $\newtilde{v}_1^{(1)} + \newtilde{v}_2^{(2)} + \newtilde{v}_3^{3} = 0$ is not met, then the mean of the divergence of the three upward-pointing triangles is in general not equal to that of the downward-pointing triangles. This leads to the creation of a so-called checkerboard pattern in the divergence. This noise also propagates from there into other fields. It is therefore important to always adhere to the condition $\newtilde{v}_1^{(1)} + \newtilde{v}_2^{(2)} + \newtilde{v}_3^{3} = 0$.
On the hexagonal grid, the reasoning can be transferred completely analogously to the relative vorticity, which leads to
\[ \begin{align} \frac{\sqrt{3}d}{4}\left(\sum_{u \in h}\zeta_u^{(t)} - \sum_{l \in h}\zeta_l^{(t)}\right) &= 6\left(\newtilde{u}_1^{(1)} + \newtilde{u}_2^{(2)} + \newtilde{u}_3^{3}\right) \end{align} \]
Failure to comply with the condition $\newtilde{u}_1^{(1)} + \newtilde{u}_2^{(2)} + \newtilde{u}_3^{3} = 0$ also leads to a checkerboard pattern on the hexagonal grid (but here in the vorticity) and must therefore also be taken into account here.
In the previous section it was shown that a checkerboard pattern exists on the hexagonal grid in the vorticity calculated on triangles if the condition $\newtilde{u}_1^{(1)} + \newtilde{u}_2^{(2)} + \newtilde{u}_3^{3} = 0$ for the wind field $\mathbf{v}$ is not met. However, in order for this condition for the wind field to be fulfilled for all time steps, this must apply to the wind field in the initial state as well as to all forcings $\mathbf{w}$. However, simulations are usually carried out with momentum diffusion to avoid a build-up of kinetic energy on the grid scale. These diffusive terms break down gradients on the grid scale and can therefore also break down an existing checkerboard pattern as long as it is not too strong. Accordingly, the requirement $\newtilde{u}_1^{(1)} + \newtilde{u}_2^{(2)} + \newtilde{u}_3^{3} = 0$ for the initial wind field $\mathbf{v}$ can be dispensed with. However, the following is important:
For all velocity tendencies $\mathbf{w}$, $\newtilde{w}_1^{(1)} + \newtilde{w}_2^{(2)} + \newtilde{w}_3^{3} = 0$ must hold.
The analogous facts apply on the triangular grid.
If one applies the Laplace operator to Eqs. (31.116) - (31.118), one obtains
\[ \begin{align} \Delta u_1 &= \frac{1}{\sqrt{3}}\left(\delta_3\Delta\newtilde{\psi}^{(2)} - \delta_2\Delta \newtilde{\psi}^{3}\right) + \delta_1\Delta\chi,\tag{31.155}\label{eq:rhombus_curl_deriv_0}\\ \Delta u_2 &= \frac{1}{\sqrt{3}}\left(\delta_1\Delta\newtilde{\psi}^{3} - \delta_3\Delta \newtilde{\psi}^{(1)}\right) + \delta_2\Delta\chi,\tag{31.156}\label{eq:rhombus_curl_deriv_1}\\ \Delta u_3 &= \frac{1}{\sqrt{3}}\left(\delta_2\Delta\newtilde{\psi}^{(1)} - \delta_1\Delta \newtilde{\psi}^{(2)}\right) + \delta_3\Delta\chi.\tag{31.157}\label{eq:rhombus_curl_deriv_2} \end{align} \]
Since discretized partial derivatives and averaging operators commute, it follows immediately
\[ \begin{align} \newtilde{\Delta u_1}^{(1)} + \newtilde{\Delta u_2}^{(2)} + \newtilde{\Delta u_3}^{(3)} = 0. \end{align} \]
If one discretizes Eqs. (31.62) - (31.64), one obtains
\[ \begin{align} \Delta u_1 &= \frac{1}{\sqrt{3}}\left(\delta_3\zeta - \delta_2\zeta\right) + \delta_1 D,\tag{31.159}\label{eq:rhombus_curl_deriv_3}\\ \Delta u_2 &= \frac{1}{\sqrt{3}}\left(\delta_1\zeta - \delta_3\zeta\right) + \delta_2 D,\tag{31.160}\label{eq:rhombus_curl_deriv_4}\\ \Delta u_3 &= \frac{1}{\sqrt{3}}\left(\delta_2\zeta - \delta_1\zeta\right) + \delta_3 D.\tag{31.161}\label{eq:rhombus_curl_deriv_5} \end{align} \]
By comparison, one finds that there are three types of vorticities in the $\mathbf{i}$ system:
\[ \begin{align} \zeta_1 \coloneqq \Delta\newtilde{\psi}^{(1)}, & {} & \zeta_2 \coloneqq \Delta\newtilde{\psi}^{(2)}, & {} & \zeta_3 \coloneqq \Delta\newtilde{\psi}^{3} \end{align} \]
It therefore holds that
\[ \begin{align} \Delta u_1 &= \frac{1}{\sqrt{3}}\left(\delta_3\zeta_2 - \delta_2\zeta_3\right) + \delta_1 D,\tag{31.163}\label{eq:rhombus_curl_deriv_6}\\ \Delta u_2 &= \frac{1}{\sqrt{3}}\left(\delta_1\zeta_3 - \delta_3\zeta_1\right) + \delta_2 D,\tag{31.164}\label{eq:rhombus_curl_deriv_7}\\ \Delta u_3 &= \frac{1}{\sqrt{3}}\left(\delta_2\zeta_1 - \delta_1\zeta_2\right) + \delta_3 D.\tag{31.165}\label{eq:rhombus_curl_deriv_8} \end{align} \]
For the divergence, according to Eq. (31.33), one has
\[ \begin{align} D &= \frac{2}{3}\left(\delta_1u_1 + \delta_2u_2 + \delta_3u_3\right). \end{align} \]
Substituting this into Eqs. (31.163) - (31.165), one obtains
\[ \begin{align} \Delta u_1 &= \frac{1}{\sqrt{3}}\left(\delta_3\zeta_2 - \delta_2\zeta_3\right) + \frac{2}{3}\left(\delta_1\delta_1u_1 + \delta_1\delta_2u_2 + \delta_1\delta_3u_3\right),\tag{31.167}\label{eq:hex_laplace_vec_0}\\ \Delta u_2 &= \frac{1}{\sqrt{3}}\left(\delta_1\zeta_3 - \delta_3\zeta_1\right) + \frac{2}{3}\left(\delta_2\delta_1u_1 + \delta_2\delta_2u_2 + \delta_2\delta_3u_3\right),\tag{31.168}\label{eq:hex_laplace_vec_1}\\ \Delta u_3 &= \frac{1}{\sqrt{3}}\left(\delta_2\zeta_1 - \delta_1\zeta_2\right) + \frac{2}{3}\left(\delta_3\delta_1u_1 + \delta_3\delta_2u_2 + \delta_3\delta_3u_3\right).\tag{31.169}\label{eq:hex_laplace_vec_2} \end{align} \]
For the vorticities $\zeta_i$, motivated by Eqs. (31.49) - (31.51), one makes the ansätze
\[ \begin{align} \zeta_1 &= \frac{2}{\sqrt{3}}\left(\delta_2u_3 - \delta_3u_2\right),\tag{31.170}\label{eq:hex_curl_rhombus_0}\\ \zeta_2 &= \frac{2}{\sqrt{3}}\left(\delta_3u_1 - \delta_1u_3\right),\tag{31.171}\label{eq:hex_curl_rhombus_1}\\ \zeta_3 &= \frac{2}{\sqrt{3}}\left(\delta_1u_2 - \delta_2u_1\right).\tag{31.172}\label{eq:hex_curl_rhombus_2} \end{align} \]
Substituting Eq. (31.170) into Eq. (31.167), one obtains
\[ \begin{align} \Delta u_1 &= \frac{2}{3}\left(\delta_3\delta_3u_1 - \delta_3\delta_1u_3 - \delta_2\delta_1u_2 + \delta_2\delta_2u_1\right) + \frac{2}{3}\left(\delta_1\delta_1u_1 + \delta_1\delta_2u_2 + \delta_1\delta_3u_3\right) = \frac{2}{3}\left(\delta_1\delta_1u_1 + \delta_2\delta_2u_1 + \delta_3\delta_3u_1\right). \end{align} \]
By cyclic permutations, the analogous statement is obtained for the other components as well. The ansätze of Eqs. (31.170) - (31.172) are therefore justified. Two neighboring triangles together form a rhombus (a parallelogram); by Eq. (31.144), such a rhombus has the area
\[ \begin{align} A_r = 2A_v = 2\frac{d^2\sqrt{3}}{4} = \frac{d^2\sqrt{3}}{2}. \end{align} \]
Thus, $\zeta_i$ corresponds to the vorticity calculated using Stokes' theorem over a rhombus whose short axis of symmetry is parallel to the $i$-axis. On the hexagonal grid, the rhombuses form the natural sets on which Stokes' theorem must be evaluated. This is a deviation from the otherwise intuitive assumption that the vorticity must be calculated via the dual grid (the triangular grid).
As usual, divergence and vorticity appear in reversed roles on the triangular grid. There, the vorticity has the unique value
\[ \begin{align} \zeta = \frac{2}{3}\left(\delta_1v_1 + \delta_2v_2 + \delta_3v_3\right) \end{align} \]
However, three types of divergences occur there, each of which is calculated using Gauss's theorem over the respective rhombuses:
\[ \begin{align} D_1 &= \frac{2}{\sqrt{3}}\left(\delta_3v_2 - \delta_2v_3\right)\\ D_2 &= \frac{2}{\sqrt{3}}\left(\delta_1v_3 - \delta_3v_1\right)\\ D_3 &= \frac{2}{\sqrt{3}}\left(\delta_2v_1 - \delta_1v_2\right) \end{align} \]
This is also underscored by the discretization of the vector Laplace operator according to Eq. (B.54), which reads
\[ \begin{align} \Delta\mathbf{v} = \nabla\left(\nabla\cdot\mathbf{v}\right) - \nabla\times\left(\nabla\times\mathbf{v}\right) \end{align} \]
Applying this to the vector component $u_1$ gives
\[ \begin{align} \Delta u_1 = \delta_1D_h - \delta_{i, \perp}\zeta_t. \end{align} \]
Here $D_h$ is the divergence on hexagons and $\zeta_t$ is the vorticity on triangles. $\delta_{i, \perp}$ denotes the central difference quotient in the $\mathbf{j}_1-$direction. For this one has
\[ \begin{align} -\delta_{i, \perp}\zeta_t &= \frac{\sqrt{3}}{d}\left(\zeta_{t, l} - \zeta_{t, u}\right). \end{align} \]
From Eq. (31.163) follows
\[ \begin{align} \Delta u_1 &= \delta_1D_h + \frac{\sqrt{3}}{d}\left(\zeta_{t, l} - \zeta_{t, u}\right) = \delta_1 D_h + \frac{1}{\sqrt{3}}\left(\delta_3\zeta_2 - \delta_2\zeta_3\right)\nonumber\\ & \Leftrightarrow \sqrt{3}\frac{\zeta_{t, l} - \zeta_{t, u}}{d} = \frac{\sqrt{3}}{3}\left(\delta_3\zeta_2 - \delta_2\zeta_3\right) \Leftrightarrow \sqrt{3}\frac{\zeta_{t, l} - \zeta_{t, u}}{d} = \frac{\sqrt{3}}{3}\frac{\zeta_{2, l} - \zeta_{2, u} - \zeta_{3, u} + \zeta_{3, l}}{d}\nonumber\\ & \Leftrightarrow \sqrt{3}\frac{\zeta_{t, l} - \zeta_{t, u}}{d} = \frac{\sqrt{3}}{3}\frac{\zeta_{2, l} + \zeta_{3, l} - \left(\zeta_{2, u} - \zeta_{3, u}\right)}{d}\nonumber\\ & \Leftrightarrow \sqrt{3}\frac{\zeta_{t, l} - \zeta_{t, u}}{d} = \frac{\sqrt{3}}{3}\frac{\zeta_{1, l} + \zeta_{2, l} + \zeta_{3, l} - \left(\zeta_{1, u} + \zeta_{2, u} - \zeta_{3, u}\right)}{d}. \end{align} \]
The last equivalence transformation follows from the fact $\zeta_{1, l} = \zeta_{1, u}$. This is satisfied by the definitions
\[ \begin{align} \zeta_{t, u} &\coloneqq \frac{\zeta_{1, u} + \zeta_{2, u} + \zeta_{3, u}}{3},\\ \zeta_{t, l} &\coloneqq \frac{\zeta_{1, l} + \zeta_{2, l} + \zeta_{3, l}}{3} \end{align} \]
The rotation on a triangle can therefore be calculated as the average of the three rhombuses that overlap with this triangle.
On the hexagonal C-grid, gradient fields $\mathbf{u} = \nabla\psi$ satisfy the condition
\[ \begin{align} \newtilde{u}_1^{(1)} + \newtilde{u}_2^{(2)} + \newtilde{u}_3^{3} = 0. \end{align} \]
On the triangular C-grid, a gradient field $\mathbf{v} = \nabla\psi$ naturally has components parallel to the basis elements $\mathbf{j}_i$. This does not satisfy
\[ \begin{align} \newtilde{v}_1^{(1)} + \newtilde{v}_2^{(2)} + \newtilde{v}_3^{3} = 0 \end{align} \]
as one can easily verify. Therefore, according to Eq. (31.153), gradient fields on the triangular C-grid produce a checkerboard pattern in the divergence, which would require unphysical numerical stabilizers such as divergence damping. This grid was therefore rightly excluded in Sect. 30.9.1. This finding was first put forward by Gassmann in [24].
The generalized Coriolis term
\[ \begin{align} \mathbf{v}\times\etabi = \mathbf{v}\times\left(\mathbf{f} + \zetabi\right) = \mathbf{v}\times\left(\mathbf{f} + \nabla\times\mathbf{v}\right) \end{align} \]
requires special attention on a C-grid. In particular, the following three questions must be answered:
Over which polygon must Stokes' theorem be evaluated to calculate the vorticity?
How must the tangential wind component on an edge be reconstructed from the normal wind components on the surrounding edges?
By what vorticity must the tangential wind component ultimately be multiplied to obtain a component of $\mathbf{v}\times\etabi$?
In this chapter, these questions will first be considered in the light of the shallow water equations (SWEs). Without friction, the SWEs read
\[ \begin{align} \frac{\partial\mathbf{v}}{\partial t} & \stackrel{\href{ch-12-important-approximations.html#eq:swe_0}{\text{Eq. (13.171)}}}{=} -g\nabla\left(h + b\right) - f\mathbf{k}\times\mathbf{v} \textcolor{red}{- \nabla k - \zeta\mathbf{k}\times\mathbf{v}} + \mu\mathbf{k},\tag{31.188}\label{eq:swe_gen_cori_0}\\ \frac{\partial h}{\partial t} & \stackrel{\href{ch-12-important-approximations.html#eq:swe_1}{\text{Eq. (13.172)}}}{=} -\nabla\cdot\left(h\mathbf{v}\right), \end{align} \]
where $\mathbf{v}$ is the velocity, $g$ is the gravitational acceleration, $h$ is the deflection of the surface from the rest position, $f$ is the Coriolis parameter, $\mathbf{k}$ is the vertical unit vector, $k \coloneqq \frac{1}{2}\mathbf{v}^2$ is the specific kinetic energy, $\zeta\coloneqq\mathbf{k}\cdot\left(\nabla\times\mathbf{v}\right)$ the relative vorticity and $\mu$ a Lagrange multiplier which, in the case of a curved set such as a sphere, ensures that the right-hand side of the equation is tangent to that set. The red terms are the momentum advection terms, which are also the nonlinear terms. Thus, the linearized momentum equation of the SWEs reads
\[ \begin{align} \frac{\partial\mathbf{v}}{\partial t} = -g\nabla\left(h + b\right) - f\mathbf{k}\times\mathbf{v}. \end{align} \]
Projecting Eq. (31.188) onto a normal vector $\mathbf{n}$ of an edge, one obtains
\[ \begin{align} \frac{\partial v_n}{\partial t} = -g\frac{\partial\left(h + b\right)}{\partial n} - f\mathbf{n}\cdot\left(\mathbf{k}\times\mathbf{v}\right) - \frac{\partial k}{\partial n} - \zeta\mathbf{n}\cdot\left(\mathbf{k}\times\mathbf{v}\right), \end{align} \]
where
\[ \begin{align} v_n \coloneqq \mathbf{n}\cdot\mathbf{v} \end{align} \]
is the normal velocity component at the edge. $\partial/\partial n$ denotes partial derivatives perpendicular to the edge. With Eq. (A.161) one obtains
\[ \begin{align} \mathbf{n}\cdot\left(\mathbf{k}\times\mathbf{v}\right) &= \mathbf{v}\cdot\left(\mathbf{n}\times\mathbf{k}\right) \stackrel{\href{ch-39-derivations-of-some-mathematical-formula.html#eq:cross_product_anticommutative}{\text{Eq. (A.151)}}}{=} -\mathbf{v}\cdot\left(\mathbf{k}\times\mathbf{n}\right) = -\mathbf{v}\cdot\mathbf{t} = -v_t,\tag{31.193}\label{eq:v_t_c-grid} \end{align} \]
where $\mathbf{t} \coloneqq \mathbf{k}\times\mathbf{n}$ is the tangential unit vector on the edge and
\[ \begin{align} v_t \coloneqq \mathbf{v}\cdot\mathbf{t} \end{align} \]
is the tangential velocity component at the edge. Thus one has
\[ \begin{align} \frac{\partial v_n}{\partial t} = -g\frac{\partial\left(h + b\right)}{\partial n} + \left(f + \zeta\right)v_t - \frac{\partial k}{\partial n}. \end{align} \]
The linear part of this is
\[ \begin{align} \frac{\partial v_n}{\partial t} = -g\frac{\partial\left(h + b\right)}{\partial n} + fv_t. \end{align} \]
The vorticity and divergence dynamics of the SWEs can be summarized in Eqs. (13.172) and (15.170):
\[ \begin{align} \frac{\partial h}{\partial t} + \nabla\cdot\left(h\mathbf{v}\right) = 0, & {} & \frac{\partial\left(hq\right)}{\partial t} + \nabla\cdot\left(hq\mathbf{v}\right) = 0 \end{align} \]
From the second equation (this is a form of the barotropic vorticity equation) one can conclude
\[ \begin{align} \frac{d}{dt}\int_{4\pi}hqd\omega &= \int_{4\pi}\frac{\partial\left(hq\right)}{\partial t}d\omega = -\int_{4\pi}\nabla\cdot\left(hq\mathbf{v}\right)d\omega \end{align} \]
The integral over the divergence of the potential vorticity flux density over the unit sphere can be evaluated by adding a third dimension in the form of a thin spherical shell of thickness $\Delta r$. No potential vorticity flows through the vertical boundaries of this set, since $w = 0$ holds in the SWEs. It thus follows
\[ \begin{align} \frac{d}{dt}\int_{4\pi}hqd\omega &= 0. \end{align} \]
The global integral of the incompressible potential vorticity $q$ is therefore constant.
The vorticity flux term of the shallow water equations reads
\[ \begin{align} -\left(f + \zeta\right)\mathbf{k}\times\mathbf{v}. \end{align} \]
In this section, a discretization is to be found that satisfies the Thuburn condition and thus, according to Eq. (31.141), also possesses a geostrophic mode. Projecting this onto the normal vector $\mathbf{n}$ of an edge, one obtains, according to Eq. (31.193),
\[ \begin{align} -\left[\left(f + \zeta\right)\mathbf{k}\times\mathbf{v}\right]\cdot\mathbf{n} &= \left(f + \zeta\right)v_t. \end{align} \]
This term is now linearized in the form
\[ \begin{align} \left(f + \zeta\right)v_t \to f_0v_t + \zeta U_t, \end{align} \]
where
\[ \begin{align} U_t \coloneqq \left(\mathbf{k}\times\mathbf{n}\right)\cdot U\mathbf{i} \end{align} \]
is the projection of a homogeneous zonal background flow $U$ onto the tangential direction. One now follows the notation from Chap. 31. There, the discretization
\[ \begin{align} f_0v_t \to \frac{f_0}{\sqrt{3}}\left(\newtilde{u}_2^{(3)} - \newtilde{u}_3^{(2)}\right) \end{align} \]
was derived for the term $fv_t$ according to Eq. (31.134). According to Eq. (31.137), this satisfies the Thuburn condition Eq. (31.102). For the term $\zeta U_t$, with
\[ \begin{align} U_{t,1} &= 0,\tag{31.205}\label{eq:vorticity_flux_c-grid_ansatz_0_0}\\ U_{t,2} &= -\zeta_2\frac{\sqrt{3}}{2}U,\tag{31.206}\label{eq:vorticity_flux_c-grid_ansatz_0_1}\\ U_{t,3} &= -\zeta_3\frac{\sqrt{3}}{2}U\tag{31.207}\label{eq:vorticity_flux_c-grid_ansatz_0_2} \end{align} \]
at the edges in the $x_i-$direction, one initially sets
\[ \begin{align} \left(\zeta U_t\right)_1 &= 0,\\ \left(\zeta U_t\right)_2 &= -\zeta_2\frac{\sqrt{3}}{2}U,\\ \left(\zeta U_t\right)_3 &= -\zeta_3\frac{\sqrt{3}}{2}U \end{align} \]
These terms must also satisfy the Thuburn condition. To check this, one computes
\[ \begin{align} \newtilde{\left(\zeta U_t\right)_1}^{(1)} + \newtilde{\left(\zeta U_t\right)_2}^{(2)} + \newtilde{\left(\zeta U_t\right)_3}^{(3)} &= -\frac{\sqrt{3}}{2}U\left(\newtilde{\zeta_2}^{(2)} + \newtilde{\zeta_3}^{(3)}\right) \stackrel{\text{in general}}{\not=} 0. \end{align} \]
Inspired by Eq. (31.134), one modifies the ansatz of Eqs. (31.205) - (31.207) to
\[ \begin{align} \left(\zeta U_t\right)_1 &= \frac{1}{\sqrt{3}}\left(-\zeta_2\frac{U}{2} + \zeta_3\frac{U}{2}\right),\\ \left(\zeta U_t\right)_2 &= \frac{1}{\sqrt{3}}\left(-\zeta_3\frac{U}{2} + \zeta_1U\right),\\ \left(\zeta U_t\right)_3 &= \frac{1}{\sqrt{3}}\left(-\zeta_1U + \zeta_2\frac{U}{2}\right). \end{align} \]
From this it follows
\[ \begin{align} \newtilde{\left(\zeta U_t\right)_1}^{(1)} + \newtilde{\left(\zeta U_t\right)_2}^{(2)} + \newtilde{\left(\zeta U_t\right)_3}^{(3)} &\propto \newtilde{\zeta_3 - \zeta_2}^{(1)} + \newtilde{2\zeta_1 - \zeta_3}^{(2)} + \newtilde{\zeta_2 - 2\zeta_1}^{(3)}\nonumber\\ = -\newtilde{\zeta_2}^{(1)} + \newtilde{\zeta_3}^{(1)} - \newtilde{\zeta_3}^{(2)} + 2\newtilde{\zeta_1}^{(2)} - 2\newtilde{\zeta_1}^{(3)} + \newtilde{\zeta_2}^{(3)} &\stackrel{\text{in general}}{\not=} 0. \end{align} \]
This ansatz also does not satisfy the Thuburn condition. Rather, a preparatory averaging of the vorticities is required, which leads to the ansatz
\[ \begin{align} \left(\zeta U_t\right)_1 &= \frac{1}{\sqrt{3}}\left(-\newoverline{\zeta_2}^{(3)}\frac{U}{2} + \newoverline{\zeta_3}^{(2)}\frac{U}{2}\right),\\ \left(\zeta U_t\right)_2 &= \frac{1}{\sqrt{3}}\left(-\newoverline{\zeta_3}^{(1)}\frac{U}{2} + \newoverline{\zeta_1}^{(3)}U\right),\\ \left(\zeta U_t\right)_3 &= \frac{1}{\sqrt{3}}\left(-\newoverline{\zeta_1}^{(2)}U + \newoverline{\zeta_2}^{(1)}\frac{U}{2}\right). \end{align} \]
This satisfies the Thuburn condition, since averaging operators commute. An operator other than the simple averaging operator $\newoverline{\zeta_i}^{(j)}$ may also be used. With $\zeta \to \eta = \zeta + f$, the complete vorticity flux term on the regular C-grid can be expressed in this way.
It is not trivial how the vorticity flux term is to be extended to the deformed grid. For this, one first starts from the f-plane in order to be able to assume a homogeneous vorticity $f_0$. Instead of attempting to generalize the Thuburn condition to the deformed grid, the closely related requirement that a geostrophic mode should exist is used as the basic condition. Furthermore, the discretization should be energy-conserving.
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{31.219}\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. (31.219), 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{31.233}\label{eq:stress_tensor_hex_regular} \end{align} \]
does not destroy any of the verified statements.
The divergence of the tensor Eq. (31.233) does not yet satisfy the Thuburn condition. Motivated by Sect. 31.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{#eq:div_3elements2d}{\text{Eq. (31.33)}}}{=} -\frac{2}{3}D + \frac{4}{3}\delta_1u_1 \end{align} \]
For the shear deformations, according to Eq. (31.233)
\[ \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{#eq:ddy_hex}{\text{Eq. (31.31)}}}{=} \frac{1}{3\sqrt{3}}\left(\delta_2 - \delta_3\right)\left(F_{1(e)} + F_{2(e)} + F_{3(e)}\right)\tag{31.236}\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. (31.236) 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{#eq:hex_curl_rhombus_1}{\text{Eq. (31.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{#eq:hex_curl_rhombus_2}{\text{Eq. (31.172)}}}{=} \zeta_3 + \frac{4}{\sqrt{3}}\delta_2u_1. \end{align} \]
Analogously to Sect. 31.6.1, one now uses the averaged shear deformations, which leads to Eq. (31.236). The procedure for calculating the gradient of the component $u_1$ is presented in Sect. 31.8.2.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.