[2282] | 1 | % ================================================================ |
---|
[3218] | 2 | % Iso-neutral diffusion : |
---|
[2282] | 3 | % ================================================================ |
---|
[2285] | 4 | \chapter{Griffies's iso-neutral diffusion} |
---|
[3218] | 5 | \label{sec:triad} |
---|
[2282] | 6 | \minitoc |
---|
| 7 | |
---|
[2285] | 8 | \section{Griffies's formulation of iso-neutral diffusion} |
---|
[3218] | 9 | \label{sec:triad:iso} |
---|
[2282] | 10 | |
---|
[3218] | 11 | We define a scheme inspired by \citet{Griffies_al_JPO98}, but formulated within the \NEMO |
---|
[2282] | 12 | framework, using scale factors rather than grid-sizes. |
---|
| 13 | |
---|
[3218] | 14 | \subsection{The iso-neutral diffusion operator} |
---|
| 15 | The iso-neutral second order tracer diffusive operator for small |
---|
| 16 | angles between iso-neutral surfaces and geopotentials is given by |
---|
| 17 | \eqref{Eq_PE_iso_tensor}: |
---|
| 18 | \begin{subequations} \label{eq:triad:PE_iso_tensor} |
---|
| 19 | \begin{equation} |
---|
| 20 | D^{lT}=-\Div\vect{f}^{lT}\equiv |
---|
| 21 | -\frac{1}{e_1e_2e_3}\left[\pd{i}\left (f_1^{lT}e_2e_3\right) + |
---|
| 22 | \pd{j}\left (f_2^{lT}e_2e_3\right) + \pd{k}\left (f_3^{lT}e_1e_2\right)\right], |
---|
| 23 | \end{equation} |
---|
| 24 | where the diffusive flux per unit area of physical space |
---|
| 25 | \begin{equation} |
---|
| 26 | \vect{f}^{lT}=-\Alt\Re\cdot\grad T, |
---|
| 27 | \end{equation} |
---|
| 28 | \begin{equation} |
---|
| 29 | \label{eq:triad:PE_iso_tensor:c} |
---|
| 30 | \mbox{with}\quad \;\;\Re = |
---|
| 31 | \begin{pmatrix} |
---|
| 32 | 1&0&-r_1\mystrut \\ |
---|
| 33 | 0&1&-r_2\mystrut \\ |
---|
| 34 | -r_1&-r_2&r_1 ^2+r_2 ^2\mystrut |
---|
| 35 | \end{pmatrix} |
---|
| 36 | \quad \text{and} \quad\grad T= |
---|
| 37 | \begin{pmatrix} |
---|
| 38 | \frac{1}{e_1}\pd[T]{i}\mystrut \\ |
---|
| 39 | \frac{1}{e_2}\pd[T]{j}\mystrut \\ |
---|
| 40 | \frac{1}{e_3}\pd[T]{k}\mystrut |
---|
| 41 | \end{pmatrix}. |
---|
| 42 | \end{equation} |
---|
| 43 | \end{subequations} |
---|
| 44 | % \left( {{\begin{array}{*{20}c} |
---|
| 45 | % 1 \hfill & 0 \hfill & {-r_1 } \hfill \\ |
---|
| 46 | % 0 \hfill & 1 \hfill & {-r_2 } \hfill \\ |
---|
| 47 | % {-r_1 } \hfill & {-r_2 } \hfill & {r_1 ^2+r_2 ^2} \hfill \\ |
---|
| 48 | % \end{array} }} \right) |
---|
| 49 | Here \eqref{Eq_PE_iso_slopes} |
---|
[2282] | 50 | \begin{align*} |
---|
| 51 | r_1 &=-\frac{e_3 }{e_1 } \left( \frac{\partial \rho }{\partial i} |
---|
| 52 | \right) |
---|
| 53 | \left( {\frac{\partial \rho }{\partial k}} \right)^{-1} \\ |
---|
| 54 | &=-\frac{e_3 }{e_1 } \left( -\alpha\frac{\partial T }{\partial i} + |
---|
| 55 | \beta\frac{\partial S }{\partial i} \right) \left( |
---|
| 56 | -\alpha\frac{\partial T }{\partial k} + \beta\frac{\partial S |
---|
| 57 | }{\partial k} \right)^{-1} |
---|
| 58 | \end{align*} |
---|
[3218] | 59 | is the $i$-component of the slope of the iso-neutral surface relative to the computational |
---|
| 60 | surface, and $r_2$ is the $j$-component. |
---|
[2282] | 61 | |
---|
[3218] | 62 | We will find it useful to consider the fluxes per unit area in $i,j,k$ |
---|
| 63 | space; we write |
---|
| 64 | \begin{equation} |
---|
| 65 | \label{eq:triad:Fijk} |
---|
| 66 | \vect{F}_{\mathrm{iso}}=\left(f_1^{lT}e_2e_3, f_2^{lT}e_1e_3, f_3^{lT}e_1e_2\right). |
---|
| 67 | \end{equation} |
---|
| 68 | Additionally, we will sometimes write the contributions towards the |
---|
| 69 | fluxes $\vect{f}$ and $\vect{F}_\mathrm{iso}$ from the component |
---|
| 70 | $R_{ij}$ of $\Re$ as $f_{ij}$, $F_{\mathrm{iso}\: ij}$, with |
---|
| 71 | $f_{ij}=R_{ij}e_i^{-1}\partial T/\partial x_i$ (no summation) etc. |
---|
| 72 | |
---|
| 73 | The off-diagonal terms of the small angle diffusion tensor |
---|
| 74 | \eqref{Eq_PE_iso_tensor}, \eqref{eq:triad:PE_iso_tensor:c} produce skew-fluxes along the |
---|
| 75 | $i$- and $j$-directions resulting from the vertical tracer gradient: |
---|
| 76 | \begin{align} |
---|
| 77 | \label{eq:triad:i13c} |
---|
| 78 | f_{13}=&+\Alt r_1\frac{1}{e_3}\frac{\partial T}{\partial k},\qquad f_{23}=+\Alt r_2\frac{1}{e_3}\frac{\partial T}{\partial k}\\ |
---|
| 79 | \intertext{and in the k-direction resulting from the lateral tracer gradients} |
---|
| 80 | \label{eq:triad:i31c} |
---|
| 81 | f_{31}+f_{32}=& \Alt r_1\frac{1}{e_1}\frac{\partial T}{\partial i}+\Alt r_2\frac{1}{e_1}\frac{\partial T}{\partial i} |
---|
| 82 | \end{align} |
---|
| 83 | |
---|
| 84 | The vertical diffusive flux associated with the $_{33}$ |
---|
[2282] | 85 | component of the small angle diffusion tensor is |
---|
| 86 | \begin{equation} |
---|
[3218] | 87 | \label{eq:triad:i33c} |
---|
| 88 | f_{33}=-\Alt(r_1^2 +r_2^2) \frac{1}{e_3}\frac{\partial T}{\partial k}. |
---|
[2282] | 89 | \end{equation} |
---|
| 90 | |
---|
| 91 | Since there are no cross terms involving $r_1$ and $r_2$ in the above, we can |
---|
[3218] | 92 | consider the iso-neutral diffusive fluxes separately in the $i$-$k$ and $j$-$k$ |
---|
[2282] | 93 | planes, just adding together the vertical components from each |
---|
[3218] | 94 | plane. The following description will describe the fluxes on the $i$-$k$ |
---|
[2282] | 95 | plane. |
---|
| 96 | |
---|
[3218] | 97 | There is no natural discretization for the $i$-component of the |
---|
| 98 | skew-flux, \eqref{eq:triad:i13c}, as |
---|
| 99 | although it must be evaluated at $u$-points, it involves vertical |
---|
[2282] | 100 | gradients (both for the tracer and the slope $r_1$), defined at |
---|
[3218] | 101 | $w$-points. Similarly, the vertical skew flux, \eqref{eq:triad:i31c}, is evaluated at |
---|
| 102 | $w$-points but involves horizontal gradients defined at $u$-points. |
---|
[2282] | 103 | |
---|
| 104 | \subsection{The standard discretization} |
---|
| 105 | The straightforward approach to discretize the lateral skew flux |
---|
[3218] | 106 | \eqref{eq:triad:i13c} from tracer cell $i,k$ to $i+1,k$, introduced in 1995 |
---|
[2282] | 107 | into OPA, \eqref{Eq_tra_ldf_iso}, is to calculate a mean vertical |
---|
[3218] | 108 | gradient at the $u$-point from the average of the four surrounding |
---|
[2282] | 109 | vertical tracer gradients, and multiply this by a mean slope at the |
---|
[3218] | 110 | $u$-point, calculated from the averaged surrounding vertical density |
---|
| 111 | gradients. The total area-integrated skew-flux (flux per unit area in |
---|
| 112 | $ijk$ space) from tracer cell $i,k$ |
---|
| 113 | to $i+1,k$, noting that the $e_{{3}_{i+1/2}^k}$ in the area |
---|
| 114 | $e{_{3}}_{i+1/2}^k{e_{2}}_{i+1/2}i^k$ at the $u$-point cancels out with |
---|
| 115 | the $1/{e_{3}}_{i+1/2}^k$ associated with the vertical tracer |
---|
[2282] | 116 | gradient, is then \eqref{Eq_tra_ldf_iso} |
---|
| 117 | \begin{equation*} |
---|
[3218] | 118 | \left(F_u^{13} \right)_{i+\hhalf}^k = \Alts_{i+\hhalf}^k |
---|
| 119 | {e_{2}}_{i+1/2}^k \overline{\overline |
---|
[2282] | 120 | r_1} ^{\,i,k}\,\overline{\overline{\delta_k T}}^{\,i,k}, |
---|
| 121 | \end{equation*} |
---|
| 122 | where |
---|
| 123 | \begin{equation*} |
---|
| 124 | \overline{\overline |
---|
[3218] | 125 | r_1} ^{\,i,k} = -\frac{{e_{3u}}_{i+1/2}^k}{{e_{1u}}_{i+1/2}^k} |
---|
| 126 | \frac{\delta_{i+1/2} [\rho]}{\overline{\overline{\delta_k \rho}}^{\,i,k}}, |
---|
[2282] | 127 | \end{equation*} |
---|
[3218] | 128 | and here and in the following we drop the $^{lT}$ superscript from |
---|
| 129 | $\Alt$ for simplicity. |
---|
[2282] | 130 | Unfortunately the resulting combination $\overline{\overline{\delta_k |
---|
| 131 | \bullet}}^{\,i,k}$ of a $k$ average and a $k$ difference %of the tracer |
---|
| 132 | reduces to $\bullet_{k+1}-\bullet_{k-1}$, so two-grid-point oscillations are |
---|
| 133 | invisible to this discretization of the iso-neutral operator. These |
---|
| 134 | \emph{computational modes} will not be damped by this operator, and |
---|
| 135 | may even possibly be amplified by it. Consequently, applying this |
---|
| 136 | operator to a tracer does not guarantee the decrease of its |
---|
| 137 | global-average variance. To correct this, we introduced a smoothing of |
---|
| 138 | the slopes of the iso-neutral surfaces (see \S\ref{LDF}). This |
---|
[3218] | 139 | technique works for $T$ and $S$ in so far as they are active tracers |
---|
[2282] | 140 | ($i.e.$ they enter the computation of density), but it does not work |
---|
| 141 | for a passive tracer. |
---|
| 142 | \subsection{Expression of the skew-flux in terms of triad slopes} |
---|
| 143 | \citep{Griffies_al_JPO98} introduce a different discretization of the |
---|
| 144 | off-diagonal terms that nicely solves the problem. |
---|
[3218] | 145 | % Instead of multiplying the mean slope calculated at the $u$-point by |
---|
| 146 | % the mean vertical gradient at the $u$-point, |
---|
[2282] | 147 | % >>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 148 | \begin{figure}[h] \begin{center} |
---|
| 149 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_triad_fluxes} |
---|
[3218] | 150 | \caption{ \label{fig:triad:ISO_triad} |
---|
[2376] | 151 | (a) Arrangement of triads $S_i$ and tracer gradients to |
---|
[3218] | 152 | give lateral tracer flux from box $i,k$ to $i+1,k$ |
---|
[2376] | 153 | (b) Triads $S'_i$ and tracer gradients to give vertical tracer flux from |
---|
| 154 | box $i,k$ to $i,k+1$.} |
---|
[3218] | 155 | \end{center} \end{figure} |
---|
[2282] | 156 | % >>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 157 | They get the skew flux from the products of the vertical gradients at |
---|
[3218] | 158 | each $w$-point surrounding the $u$-point with the corresponding `triad' |
---|
| 159 | slope calculated from the lateral density gradient across the $u$-point |
---|
| 160 | divided by the vertical density gradient at the same $w$-point as the |
---|
| 161 | tracer gradient. See Fig.~\ref{fig:triad:ISO_triad}a, where the thick lines |
---|
[2282] | 162 | denote the tracer gradients, and the thin lines the corresponding |
---|
| 163 | triads, with slopes $s_1, \dotsc s_4$. The total area-integrated |
---|
| 164 | skew-flux from tracer cell $i,k$ to $i+1,k$ |
---|
| 165 | \begin{multline} |
---|
[3218] | 166 | \label{eq:triad:i13} |
---|
| 167 | \left( F_u^{13} \right)_{i+\frac{1}{2}}^k = \Alts_{i+1}^k a_1 s_1 |
---|
[2282] | 168 | \delta _{k+\frac{1}{2}} \left[ T^{i+1} |
---|
[3218] | 169 | \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}} + \Alts _i^k a_2 s_2 \delta |
---|
[2282] | 170 | _{k+\frac{1}{2}} \left[ T^i |
---|
| 171 | \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}} \\ |
---|
[3218] | 172 | +\Alts _{i+1}^k a_3 s_3 \delta _{k-\frac{1}{2}} \left[ T^{i+1} |
---|
| 173 | \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}} +\Alts _i^k a_4 s_4 \delta |
---|
[2282] | 174 | _{k-\frac{1}{2}} \left[ T^i \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}}, |
---|
| 175 | \end{multline} |
---|
| 176 | where the contributions of the triad fluxes are weighted by areas |
---|
[3218] | 177 | $a_1, \dotsc a_4$, and $\Alts$ is now defined at the tracer points |
---|
| 178 | rather than the $u$-points. This discretization gives a much closer |
---|
[2282] | 179 | stencil, and disallows the two-point computational modes. |
---|
| 180 | |
---|
[3218] | 181 | The vertical skew flux \eqref{eq:triad:i31c} from tracer cell $i,k$ to $i,k+1$ at the |
---|
| 182 | $w$-point $i,k+\hhalf$ is constructed similarly (Fig.~\ref{fig:triad:ISO_triad}b) |
---|
[2282] | 183 | by multiplying lateral tracer gradients from each of the four |
---|
[3218] | 184 | surrounding $u$-points by the appropriate triad slope: |
---|
[2282] | 185 | \begin{multline} |
---|
[3218] | 186 | \label{eq:triad:i31} |
---|
| 187 | \left( F_w^{31} \right) _i ^{k+\frac{1}{2}} = \Alts_i^{k+1} a_{1}' |
---|
[2282] | 188 | s_{1}' \delta _{i-\frac{1}{2}} \left[ T^{k+1} \right]/{e_{3u}}_{i-\frac{1}{2}}^{k+1} |
---|
[3218] | 189 | +\Alts_i^{k+1} a_{2}' s_{2}' \delta _{i+\frac{1}{2}} \left[ T^{k+1} \right]/{e_{3u}}_{i+\frac{1}{2}}^{k+1}\\ |
---|
| 190 | + \Alts_i^k a_{3}' s_{3}' \delta _{i-\frac{1}{2}} \left[ T^k\right]/{e_{3u}}_{i-\frac{1}{2}}^k |
---|
| 191 | +\Alts_i^k a_{4}' s_{4}' \delta _{i+\frac{1}{2}} \left[ T^k \right]/{e_{3u}}_{i+\frac{1}{2}}^k. |
---|
[2282] | 192 | \end{multline} |
---|
[3218] | 193 | |
---|
| 194 | We notate the triad slopes $s_i$ and $s'_i$ in terms of the `anchor point' $i,k$ |
---|
| 195 | (appearing in both the vertical and lateral gradient), and the $u$- and |
---|
| 196 | $w$-points $(i+i_p,k)$, $(i,k+k_p)$ at the centres of the `arms' of the |
---|
| 197 | triad as follows (see also Fig.~\ref{fig:triad:ISO_triad}): |
---|
[2282] | 198 | \begin{equation} |
---|
[3218] | 199 | \label{eq:triad:R} |
---|
| 200 | _i^k \mathbb{R}_{i_p}^{k_p} |
---|
| 201 | =-\frac{ {e_{3w}}_{\,i}^{\,k+k_p}} { {e_{1u}}_{\,i+i_p}^{\,k}} |
---|
[2282] | 202 | \ |
---|
[3218] | 203 | \frac |
---|
[2282] | 204 | {\left(\alpha / \beta \right)_i^k \ \delta_{i + i_p}[T^k] - \delta_{i + i_p}[S^k] } |
---|
| 205 | {\left(\alpha / \beta \right)_i^k \ \delta_{k+k_p}[T^i ] - \delta_{k+k_p}[S^i ] }. |
---|
| 206 | \end{equation} |
---|
| 207 | In calculating the slopes of the local neutral |
---|
| 208 | surfaces, the expansion coefficients $\alpha$ and $\beta$ are |
---|
[3218] | 209 | evaluated at the anchor points of the triad \footnote{Note that in \eqref{eq:triad:R} we use the ratio $\alpha / \beta$ |
---|
[2282] | 210 | instead of multiplying the temperature derivative by $\alpha$ and the |
---|
| 211 | salinity derivative by $\beta$. This is more efficient as the ratio |
---|
| 212 | $\alpha / \beta$ can to be evaluated directly}, while the metrics are |
---|
[3218] | 213 | calculated at the $u$- and $w$-points on the arms. |
---|
[2282] | 214 | |
---|
| 215 | % >>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 216 | \begin{figure}[h] \begin{center} |
---|
| 217 | \includegraphics[width=0.60\textwidth]{./TexFiles/Figures/Fig_qcells} |
---|
[3218] | 218 | \caption{ \label{fig:triad:qcells} |
---|
| 219 | Triad notation for quarter cells.$T$-cells are inside |
---|
| 220 | boxes, while the $i+\half,k$ $u$-cell is shaded in green and the |
---|
| 221 | $i,k+\half$ $w$-cell is shaded in pink.} |
---|
[2282] | 222 | \end{center} \end{figure} |
---|
| 223 | % >>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 224 | |
---|
[3218] | 225 | Each triad $\{_i^k\:_{i_p}^{k_p}\}$ is associated (Fig.~\ref{fig:triad:qcells}) with the quarter |
---|
| 226 | cell that is the intersection of the $i,k$ $T$-cell, the $i+i_p,k$ |
---|
| 227 | $u$-cell and the $i,k+k_p$ $w$-cell. Expressing the slopes $s_i$ and |
---|
| 228 | $s'_i$ in \eqref{eq:triad:i13} and \eqref{eq:triad:i31} in this notation, we have |
---|
[2282] | 229 | e.g.\ $s_1=s'_1={\:}_i^k \mathbb{R}_{1/2}^{1/2}$. Each triad slope $_i^k |
---|
| 230 | \mathbb{R}_{i_p}^{k_p}$ is used once (as an $s$) to calculate the |
---|
[3218] | 231 | lateral flux along its $u$-arm, at $(i+i_p,k)$, and then again as an |
---|
| 232 | $s'$ to calculate the vertical flux along its $w$-arm at |
---|
[2282] | 233 | $(i,k+k_p)$. Each vertical area $a_i$ used to calculate the lateral |
---|
| 234 | flux and horizontal area $a'_i$ used to calculate the vertical flux |
---|
[3218] | 235 | can also be identified as the area across the $u$- and $w$-arms of a |
---|
| 236 | unique triad, and we notate these areas, similarly to the triad |
---|
[2282] | 237 | slopes, as $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$, |
---|
[3218] | 238 | $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$, where e.g. in \eqref{eq:triad:i13} |
---|
| 239 | $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$, and in \eqref{eq:triad:i31} |
---|
[2282] | 240 | $a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$. |
---|
| 241 | |
---|
| 242 | \subsection{The full triad fluxes} |
---|
[3218] | 243 | A key property of iso-neutral diffusion is that it should not affect |
---|
[2282] | 244 | the (locally referenced) density. In particular there should be no |
---|
| 245 | lateral or vertical density flux. The lateral density flux disappears so long as the |
---|
| 246 | area-integrated lateral diffusive flux from tracer cell $i,k$ to |
---|
| 247 | $i+1,k$ coming from the $_{11}$ term of the diffusion tensor takes the |
---|
| 248 | form |
---|
| 249 | \begin{equation} |
---|
[3218] | 250 | \label{eq:triad:i11} |
---|
[2282] | 251 | \left( F_u^{11} \right) _{i+\frac{1}{2}} ^{k} = |
---|
[3218] | 252 | - \left( \Alts_i^{k+1} a_{1} + \Alts_i^{k+1} a_{2} + \Alts_i^k |
---|
| 253 | a_{3} + \Alts_i^k a_{4} \right) |
---|
[2282] | 254 | \frac{\delta _{i+1/2} \left[ T^k\right]}{{e_{1u}}_{\,i+1/2}^{\,k}}, |
---|
| 255 | \end{equation} |
---|
[3218] | 256 | where the areas $a_i$ are as in \eqref{eq:triad:i13}. In this case, |
---|
| 257 | separating the total lateral flux, the sum of \eqref{eq:triad:i13} and |
---|
| 258 | \eqref{eq:triad:i11}, into triad components, a lateral tracer |
---|
[2282] | 259 | flux |
---|
| 260 | \begin{equation} |
---|
[3218] | 261 | \label{eq:triad:latflux-triad} |
---|
| 262 | _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) = - \Alts_i^k{ \:}_i^k{\mathbb{A}_u}_{i_p}^{k_p} |
---|
[2282] | 263 | \left( |
---|
| 264 | \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } |
---|
| 265 | -\ {_i^k\mathbb{R}_{i_p}^{k_p}} \ |
---|
| 266 | \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } |
---|
| 267 | \right) |
---|
| 268 | \end{equation} |
---|
| 269 | can be identified with each triad. Then, because the |
---|
| 270 | same metric factors ${e_{3w}}_{\,i}^{\,k+k_p}$ and |
---|
| 271 | ${e_{1u}}_{\,i+i_p}^{\,k}$ are employed for both the density gradients |
---|
| 272 | in $ _i^k \mathbb{R}_{i_p}^{k_p}$ and the tracer gradients, the lateral |
---|
| 273 | density flux associated with each triad separately disappears. |
---|
| 274 | \begin{equation} |
---|
[3218] | 275 | \label{eq:triad:latflux-rho} |
---|
[2282] | 276 | {\mathbb{F}_u}_{i_p}^{k_p} (\rho)=-\alpha _i^k {\:}_i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) + \beta_i^k {\:}_i^k {\mathbb{F}_u}_{i_p}^{k_p} (S)=0 |
---|
| 277 | \end{equation} |
---|
| 278 | Thus the total flux $\left( F_u^{31} \right) ^i _{i,k+\frac{1}{2}} + |
---|
| 279 | \left( F_u^{11} \right) ^i _{i,k+\frac{1}{2}}$ from tracer cell $i,k$ |
---|
| 280 | to $i+1,k$ must also vanish since it is a sum of four such triad fluxes. |
---|
| 281 | |
---|
[3218] | 282 | The squared slope $r_1^2$ in the expression \eqref{eq:triad:i33c} for the |
---|
[2282] | 283 | $_{33}$ component is also expressed in terms of area-weighted |
---|
| 284 | squared triad slopes, so the area-integrated vertical flux from tracer |
---|
| 285 | cell $i,k$ to $i,k+1$ resulting from the $r_1^2$ term is |
---|
| 286 | \begin{equation} |
---|
[3218] | 287 | \label{eq:triad:i33} |
---|
[2282] | 288 | \left( F_w^{33} \right) _i^{k+\frac{1}{2}} = |
---|
[3218] | 289 | - \left( \Alts_i^{k+1} a_{1}' s_{1}'^2 |
---|
| 290 | + \Alts_i^{k+1} a_{2}' s_{2}'^2 |
---|
| 291 | + \Alts_i^k a_{3}' s_{3}'^2 |
---|
| 292 | + \Alts_i^k a_{4}' s_{4}'^2 \right)\delta_{k+\frac{1}{2}} \left[ T^{i+1} \right], |
---|
[2282] | 293 | \end{equation} |
---|
| 294 | where the areas $a'$ and slopes $s'$ are the same as in |
---|
[3218] | 295 | \eqref{eq:triad:i31}. |
---|
| 296 | Then, separating the total vertical flux, the sum of \eqref{eq:triad:i31} and |
---|
| 297 | \eqref{eq:triad:i33}, into triad components, a vertical flux |
---|
[2282] | 298 | \begin{align} |
---|
[3218] | 299 | \label{eq:triad:vertflux-triad} |
---|
[2282] | 300 | _i^k {\mathbb{F}_w}_{i_p}^{k_p} (T) |
---|
[3218] | 301 | &= \Alts_i^k{\: }_i^k{\mathbb{A}_w}_{i_p}^{k_p} |
---|
[2282] | 302 | \left( |
---|
| 303 | {_i^k\mathbb{R}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } |
---|
| 304 | -\ \left({_i^k\mathbb{R}_{i_p}^{k_p}}\right)^2 \ |
---|
| 305 | \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } |
---|
| 306 | \right) \\ |
---|
| 307 | &= - \left(\left.{ }_i^k{\mathbb{A}_w}_{i_p}^{k_p}\right/{ }_i^k{\mathbb{A}_u}_{i_p}^{k_p}\right) |
---|
[3218] | 308 | {_i^k\mathbb{R}_{i_p}^{k_p}}{\: }_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \label{eq:triad:vertflux-triad2} |
---|
[2282] | 309 | \end{align} |
---|
| 310 | may be associated with each triad. Each vertical density flux $_i^k {\mathbb{F}_w}_{i_p}^{k_p} (\rho)$ |
---|
| 311 | associated with a triad then separately disappears (because the |
---|
| 312 | lateral flux $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (\rho)$ |
---|
| 313 | disappears). Consequently the total vertical density flux $\left( F_w^{31} \right)_i ^{k+\frac{1}{2}} + |
---|
| 314 | \left( F_w^{33} \right)_i^{k+\frac{1}{2}}$ from tracer cell $i,k$ |
---|
| 315 | to $i,k+1$ must also vanish since it is a sum of four such triad |
---|
| 316 | fluxes. |
---|
| 317 | |
---|
[3218] | 318 | We can explicitly identify (Fig.~\ref{fig:triad:qcells}) the triads associated with the $s_i$, $a_i$, and $s'_i$, $a'_i$ used in the definition of |
---|
| 319 | the $u$-fluxes and $w$-fluxes in |
---|
| 320 | \eqref{eq:triad:i31}, \eqref{eq:triad:i13}, \eqref{eq:triad:i11} \eqref{eq:triad:i33} and |
---|
| 321 | Fig.~\ref{fig:triad:ISO_triad} to write out the iso-neutral fluxes at $u$- and |
---|
[2282] | 322 | $w$-points as sums of the triad fluxes that cross the $u$- and $w$-faces: |
---|
| 323 | %(Fig.~\ref{Fig_ISO_triad}): |
---|
[3218] | 324 | \begin{flalign} \label{Eq_iso_flux} \vect{F}_\mathrm{iso}(T) &\equiv |
---|
[2282] | 325 | \sum_{\substack{i_p,\,k_p}} |
---|
| 326 | \begin{pmatrix} |
---|
| 327 | {_{i+1/2-i_p}^k {\mathbb{F}_u}_{i_p}^{k_p} } (T) \\ |
---|
| 328 | \\ |
---|
| 329 | {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p} } (T) \\ |
---|
| 330 | \end{pmatrix}. |
---|
| 331 | \end{flalign} |
---|
[3218] | 332 | \subsection{Ensuring the scheme does not increase tracer variance} |
---|
| 333 | \label{sec:triad:variance} |
---|
[2282] | 334 | |
---|
[3218] | 335 | We now require that this operator should not increase the |
---|
[2282] | 336 | globally-integrated tracer variance. |
---|
| 337 | %This changes according to |
---|
| 338 | % \begin{align*} |
---|
| 339 | % &\int_D D_l^T \; T \;dv \equiv \sum_{i,k} \left\{ T \ D_l^T \ b_T \right\} \\ |
---|
[3218] | 340 | % &\equiv + \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{ |
---|
| 341 | % \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right] |
---|
[2282] | 342 | % + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right] \ T \right\} \\ |
---|
[3218] | 343 | % &\equiv - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{ |
---|
[2282] | 344 | % {_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \ \delta_{i+1/2} [T] |
---|
| 345 | % + {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \ \delta_{k+1/2} [T] \right\} \\ |
---|
| 346 | % \end{align*} |
---|
| 347 | Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ drives a lateral flux |
---|
[3218] | 348 | $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ across the $u$-point $i+i_p,k$ and |
---|
[2282] | 349 | a vertical flux $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ across the |
---|
[3218] | 350 | $w$-point $i,k+k_p$. The lateral flux drives a net rate of change of |
---|
| 351 | variance, summed over the two $T$-points $i+i_p-\half,k$ and $i+i_p+\half,k$, of |
---|
[2282] | 352 | \begin{multline} |
---|
| 353 | {b_T}_{i+i_p-1/2}^k\left(\frac{\partial T}{\partial t}T\right)_{i+i_p-1/2}^k+ |
---|
| 354 | \quad {b_T}_{i+i_p+1/2}^k\left(\frac{\partial T}{\partial |
---|
| 355 | t}T\right)_{i+i_p+1/2}^k \\ |
---|
| 356 | \begin{split} |
---|
| 357 | &= -T_{i+i_p-1/2}^k{\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \quad + \quad T_{i+i_p+1/2}^k |
---|
| 358 | {\;}_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \\ |
---|
[3218] | 359 | &={\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], \label{eq:triad:dvar_iso_i} |
---|
[2282] | 360 | \end{split} |
---|
| 361 | \end{multline} |
---|
[3218] | 362 | while the vertical flux similarly drives a net rate of change of |
---|
| 363 | variance summed over the $T$-points $i,k+k_p-\half$ (above) and |
---|
| 364 | $i,k+k_p+\half$ (below) of |
---|
[2282] | 365 | \begin{equation} |
---|
[3218] | 366 | \label{eq:triad:dvar_iso_k} |
---|
[2282] | 367 | _i^k{\mathbb{F}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i]. |
---|
| 368 | \end{equation} |
---|
| 369 | The total variance tendency driven by the triad is the sum of these |
---|
| 370 | two. Expanding $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ and |
---|
[3218] | 371 | $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ with \eqref{eq:triad:latflux-triad} and |
---|
| 372 | \eqref{eq:triad:vertflux-triad}, it is |
---|
[2282] | 373 | \begin{multline*} |
---|
[3218] | 374 | -\Alts_i^k\left \{ |
---|
[2282] | 375 | { } _i^k{\mathbb{A}_u}_{i_p}^{k_p} |
---|
| 376 | \left( |
---|
| 377 | \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } |
---|
| 378 | - {_i^k\mathbb{R}_{i_p}^{k_p}} \ |
---|
| 379 | \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }\right)\,\delta_{i+ i_p}[T^k] \right.\\ |
---|
| 380 | - \left. { } _i^k{\mathbb{A}_w}_{i_p}^{k_p} |
---|
| 381 | \left( |
---|
| 382 | \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } |
---|
| 383 | -{\:}_i^k\mathbb{R}_{i_p}^{k_p} |
---|
| 384 | \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } |
---|
| 385 | \right) {\,}_i^k\mathbb{R}_{i_p}^{k_p}\delta_{k+ k_p}[T^i] |
---|
| 386 | \right \}. |
---|
| 387 | \end{multline*} |
---|
| 388 | The key point is then that if we require |
---|
| 389 | $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$ and $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$ |
---|
| 390 | to be related to a triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$ by |
---|
| 391 | \begin{equation} |
---|
[3218] | 392 | \label{eq:triad:V-A} |
---|
[2282] | 393 | _i^k\mathbb{V}_{i_p}^{k_p} |
---|
| 394 | ={\;}_i^k{\mathbb{A}_u}_{i_p}^{k_p}\,{e_{1u}}_{\,i + i_p}^{\,k} |
---|
| 395 | ={\;}_i^k{\mathbb{A}_w}_{i_p}^{k_p}\,{e_{3w}}_{\,i}^{\,k + k_p}, |
---|
| 396 | \end{equation} |
---|
| 397 | the variance tendency reduces to the perfect square |
---|
| 398 | \begin{equation} |
---|
[3218] | 399 | \label{eq:triad:perfect-square} |
---|
| 400 | -\Alts_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p} |
---|
[2282] | 401 | \left( |
---|
| 402 | \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } |
---|
| 403 | -{\:}_i^k\mathbb{R}_{i_p}^{k_p} |
---|
| 404 | \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } |
---|
| 405 | \right)^2\leq 0. |
---|
| 406 | \end{equation} |
---|
[3218] | 407 | Thus, the constraint \eqref{eq:triad:V-A} ensures that the fluxes (\ref{eq:triad:latflux-triad}, \ref{eq:triad:vertflux-triad}) associated |
---|
[2282] | 408 | with a given slope triad $_i^k\mathbb{R}_{i_p}^{k_p}$ do not increase |
---|
| 409 | the net variance. Since the total fluxes are sums of such fluxes from |
---|
| 410 | the various triads, this constraint, applied to all triads, is |
---|
| 411 | sufficient to ensure that the globally integrated variance does not |
---|
| 412 | increase. |
---|
| 413 | |
---|
[3218] | 414 | The expression \eqref{eq:triad:V-A} can be interpreted as a discretization |
---|
[2282] | 415 | of the global integral |
---|
| 416 | \begin{equation} |
---|
[3218] | 417 | \label{eq:triad:cts-var} |
---|
| 418 | \frac{\partial}{\partial t}\int\!\half T^2\, dV = |
---|
| 419 | \int\!\mathbf{F}\cdot\nabla T\, dV, |
---|
[2282] | 420 | \end{equation} |
---|
| 421 | where, within each triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$, the |
---|
| 422 | lateral and vertical fluxes/unit area |
---|
| 423 | \[ |
---|
| 424 | \mathbf{F}=\left( |
---|
[3218] | 425 | \left.{}_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\right/{}_i^k{\mathbb{A}_u}_{i_p}^{k_p}, |
---|
| 426 | \left.{\:}_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)\right/{}_i^k{\mathbb{A}_w}_{i_p}^{k_p} |
---|
[2282] | 427 | \right) |
---|
| 428 | \] |
---|
| 429 | and the gradient |
---|
| 430 | \[\nabla T = \left( |
---|
[3218] | 431 | \left.\delta_{i+ i_p}[T^k] \right/ {e_{1u}}_{\,i + i_p}^{\,k}, |
---|
| 432 | \left.\delta_{k+ k_p}[T^i] \right/ {e_{3w}}_{\,i}^{\,k + k_p} |
---|
[2282] | 433 | \right) |
---|
| 434 | \] |
---|
| 435 | \subsection{Triad volumes in Griffes's scheme and in \NEMO} |
---|
| 436 | To complete the discretization we now need only specify the triad |
---|
| 437 | volumes $_i^k\mathbb{V}_{i_p}^{k_p}$. \citet{Griffies_al_JPO98} identify |
---|
| 438 | these $_i^k\mathbb{V}_{i_p}^{k_p}$ as the volumes of the quarter |
---|
[3218] | 439 | cells, defined in terms of the distances between $T$, $u$,$f$ and |
---|
| 440 | $w$-points. This is the natural discretization of |
---|
| 441 | \eqref{eq:triad:cts-var}. The \NEMO model, however, operates with scale |
---|
[2282] | 442 | factors instead of grid sizes, and scale factors for the quarter |
---|
| 443 | cells are not defined. Instead, therefore we simply choose |
---|
| 444 | \begin{equation} |
---|
[3218] | 445 | \label{eq:triad:V-NEMO} |
---|
[2282] | 446 | _i^k\mathbb{V}_{i_p}^{k_p}=\quarter {b_u}_{i+i_p}^k, |
---|
| 447 | \end{equation} |
---|
[3218] | 448 | as a quarter of the volume of the $u$-cell inside which the triad |
---|
[2282] | 449 | quarter-cell lies. This has the nice property that when the slopes |
---|
| 450 | $\mathbb{R}$ vanish, the lateral flux from tracer cell $i,k$ to |
---|
| 451 | $i+1,k$ reduces to the classical form |
---|
| 452 | \begin{equation} |
---|
[3218] | 453 | \label{eq:triad:lat-normal} |
---|
| 454 | -\overline\Alts_{\,i+1/2}^k\; |
---|
[2282] | 455 | \frac{{b_u}_{i+1/2}^k}{{e_{1u}}_{\,i + i_p}^{\,k}} |
---|
| 456 | \;\frac{\delta_{i+ 1/2}[T^k] }{{e_{1u}}_{\,i + i_p}^{\,k}} |
---|
[3218] | 457 | = -\overline\Alts_{\,i+1/2}^k\;\frac{{e_{1w}}_{\,i + 1/2}^{\,k}\:{e_{1v}}_{\,i + 1/2}^{\,k}\;\delta_{i+ 1/2}[T^k]}{{e_{1u}}_{\,i + 1/2}^{\,k}}. |
---|
[2282] | 458 | \end{equation} |
---|
[3218] | 459 | In fact if the diffusive coefficient is defined at $u$-points, so that |
---|
| 460 | we employ $\Alts_{i+i_p}^k$ instead of $\Alts_i^k$ in the definitions of the |
---|
| 461 | triad fluxes \eqref{eq:triad:latflux-triad} and \eqref{eq:triad:vertflux-triad}, |
---|
[2282] | 462 | we can replace $\overline{A}_{\,i+1/2}^k$ by $A_{i+1/2}^k$ in the above. |
---|
| 463 | |
---|
| 464 | \subsection{Summary of the scheme} |
---|
[3218] | 465 | The iso-neutral fluxes at $u$- and |
---|
| 466 | $w$-points are the sums of the triad fluxes that cross the $u$- and |
---|
| 467 | $w$-faces \eqref{Eq_iso_flux}: |
---|
| 468 | \begin{subequations}\label{eq:triad:alltriadflux} |
---|
| 469 | \begin{flalign}\label{eq:triad:vect_isoflux} |
---|
| 470 | \vect{F}_\mathrm{iso}(T) &\equiv |
---|
| 471 | \sum_{\substack{i_p,\,k_p}} |
---|
| 472 | \begin{pmatrix} |
---|
| 473 | {_{i+1/2-i_p}^k {\mathbb{F}_u}_{i_p}^{k_p} } (T) \\ |
---|
| 474 | \\ |
---|
| 475 | {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p} } (T) |
---|
| 476 | \end{pmatrix}, |
---|
| 477 | \end{flalign} |
---|
| 478 | where \eqref{eq:triad:latflux-triad}: |
---|
| 479 | \begin{align} |
---|
| 480 | \label{eq:triad:triadfluxu} |
---|
| 481 | _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) &= - \Alts_i^k{ |
---|
| 482 | \:}\frac{{{}_i^k\mathbb{V}}_{i_p}^{k_p}}{{e_{1u}}_{\,i + i_p}^{\,k}} |
---|
| 483 | \left( |
---|
| 484 | \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } |
---|
| 485 | -\ {_i^k\mathbb{R}_{i_p}^{k_p}} \ |
---|
| 486 | \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } |
---|
| 487 | \right),\\ |
---|
| 488 | \intertext{and} |
---|
| 489 | _i^k {\mathbb{F}_w}_{i_p}^{k_p} (T) |
---|
| 490 | &= \Alts_i^k{\: }\frac{{{}_i^k\mathbb{V}}_{i_p}^{k_p}}{{e_{3w}}_{\,i}^{\,k+k_p}} |
---|
| 491 | \left( |
---|
| 492 | {_i^k\mathbb{R}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } |
---|
| 493 | -\ \left({_i^k\mathbb{R}_{i_p}^{k_p}}\right)^2 \ |
---|
| 494 | \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } |
---|
| 495 | \right),\label{eq:triad:triadfluxw} |
---|
| 496 | \end{align} |
---|
| 497 | with \eqref{eq:triad:V-NEMO} |
---|
| 498 | \begin{equation} |
---|
| 499 | \label{eq:triad:V-NEMO2} |
---|
| 500 | _i^k{\mathbb{V}}_{i_p}^{k_p}=\quarter {b_u}_{i+i_p}^k. |
---|
| 501 | \end{equation} |
---|
| 502 | \end{subequations} |
---|
| 503 | |
---|
| 504 | The divergence of the expression \eqref{Eq_iso_flux} for the fluxes gives the iso-neutral diffusion tendency at |
---|
[2282] | 505 | each tracer point: |
---|
[3218] | 506 | \begin{equation} \label{eq:triad:iso_operator} D_l^T = \frac{1}{b_T} |
---|
[2282] | 507 | \sum_{\substack{i_p,\,k_p}} \left\{ \delta_{i} \left[{_{i+1/2-i_p}^k |
---|
| 508 | {\mathbb{F}_u }_{i_p}^{k_p}} \right] + \delta_{k} \left[ |
---|
| 509 | {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right] \right\} |
---|
| 510 | \end{equation} |
---|
| 511 | where $b_T= e_{1T}\,e_{2T}\,e_{3T}$ is the volume of $T$-cells. |
---|
| 512 | The diffusion scheme satisfies the following six properties: |
---|
| 513 | \begin{description} |
---|
| 514 | \item[$\bullet$ horizontal diffusion] The discretization of the |
---|
[3218] | 515 | diffusion operator recovers \eqref{eq:triad:lat-normal} the traditional five-point Laplacian in |
---|
[2282] | 516 | the limit of flat iso-neutral direction : |
---|
[3218] | 517 | \begin{equation} \label{eq:triad:iso_property0} D_l^T = \frac{1}{b_T} \ |
---|
[2282] | 518 | \delta_{i} \left[ \frac{e_{2u}\,e_{3u}}{e_{1u}} \; |
---|
[3218] | 519 | \overline\Alts^{\,i} \; \delta_{i+1/2}[T] \right] \qquad |
---|
[2282] | 520 | \text{when} \quad { _i^k \mathbb{R}_{i_p}^{k_p} }=0 |
---|
| 521 | \end{equation} |
---|
| 522 | |
---|
| 523 | \item[$\bullet$ implicit treatment in the vertical] Only tracer values |
---|
| 524 | associated with a single water column appear in the expression |
---|
[3218] | 525 | \eqref{eq:triad:i33} for the $_{33}$ fluxes, vertical fluxes driven by |
---|
[2282] | 526 | vertical gradients. This is of paramount importance since it means |
---|
[3218] | 527 | that a time-implicit algorithm can be used to solve the vertical |
---|
| 528 | diffusion equation. This is necessary |
---|
| 529 | since the vertical eddy |
---|
[2282] | 530 | diffusivity associated with this term, |
---|
| 531 | \begin{equation} |
---|
[3218] | 532 | \frac{1}{b_w}\sum_{\substack{i_p, \,k_p}} \left\{ |
---|
| 533 | {\:}_i^k\mathbb{V}_{i_p}^{k_p} \: \Alts_i^k \: \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2 |
---|
| 534 | \right\} = |
---|
| 535 | \frac{1}{4b_w}\sum_{\substack{i_p, \,k_p}} \left\{ |
---|
| 536 | {b_u}_{i+i_p}^k\: \Alts_i^k \: \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2 |
---|
[2282] | 537 | \right\}, |
---|
| 538 | \end{equation} |
---|
| 539 | (where $b_w= e_{1w}\,e_{2w}\,e_{3w}$ is the volume of $w$-cells) can be quite large. |
---|
| 540 | |
---|
| 541 | \item[$\bullet$ pure iso-neutral operator] The iso-neutral flux of |
---|
| 542 | locally referenced potential density is zero. See |
---|
[3218] | 543 | \eqref{eq:triad:latflux-rho} and \eqref{eq:triad:vertflux-triad2}. |
---|
[2282] | 544 | |
---|
| 545 | \item[$\bullet$ conservation of tracer] The iso-neutral diffusion |
---|
| 546 | conserves tracer content, $i.e.$ |
---|
[3218] | 547 | \begin{equation} \label{eq:triad:iso_property1} \sum_{i,j,k} \left\{ D_l^T \ |
---|
[2282] | 548 | b_T \right\} = 0 |
---|
| 549 | \end{equation} |
---|
| 550 | This property is trivially satisfied since the iso-neutral diffusive |
---|
| 551 | operator is written in flux form. |
---|
| 552 | |
---|
| 553 | \item[$\bullet$ no increase of tracer variance] The iso-neutral diffusion |
---|
| 554 | does not increase the tracer variance, $i.e.$ |
---|
[3218] | 555 | \begin{equation} \label{eq:triad:iso_property2} \sum_{i,j,k} \left\{ T \ D_l^T |
---|
[2282] | 556 | \ b_T \right\} \leq 0 |
---|
| 557 | \end{equation} |
---|
| 558 | The property is demonstrated in |
---|
[3218] | 559 | \S\ref{sec:triad:variance} above. It is a key property for a diffusion |
---|
| 560 | term. It means that it is also a dissipation term, $i.e.$ it |
---|
| 561 | dissipates the square of the quantity on which it is applied. It |
---|
[2282] | 562 | therefore ensures that, when the diffusivity coefficient is large |
---|
[3218] | 563 | enough, the field on which it is applied becomes free of grid-point |
---|
[2282] | 564 | noise. |
---|
| 565 | |
---|
| 566 | \item[$\bullet$ self-adjoint operator] The iso-neutral diffusion |
---|
| 567 | operator is self-adjoint, $i.e.$ |
---|
[3218] | 568 | \begin{equation} \label{eq:triad:iso_property3} \sum_{i,j,k} \left\{ S \ D_l^T |
---|
[2282] | 569 | \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\} |
---|
| 570 | \end{equation} |
---|
| 571 | In other word, there is no need to develop a specific routine from |
---|
| 572 | the adjoint of this operator. We just have to apply the same |
---|
| 573 | routine. This property can be demonstrated similarly to the proof of |
---|
| 574 | the `no increase of tracer variance' property. The contribution by a |
---|
[3218] | 575 | single triad towards the left hand side of \eqref{eq:triad:iso_property3}, can |
---|
| 576 | be found by replacing $\delta[T]$ by $\delta[S]$ in \eqref{eq:triad:dvar_iso_i} |
---|
| 577 | and \eqref{eq:triad:dvar_iso_k}. This results in a term similar to |
---|
| 578 | \eqref{eq:triad:perfect-square}, |
---|
[2282] | 579 | \begin{equation} |
---|
[3218] | 580 | \label{eq:triad:TScovar} |
---|
| 581 | - \Alts_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p} |
---|
[2282] | 582 | \left( |
---|
| 583 | \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } |
---|
| 584 | -{\:}_i^k\mathbb{R}_{i_p}^{k_p} |
---|
| 585 | \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } |
---|
| 586 | \right) |
---|
| 587 | \left( |
---|
| 588 | \frac{ \delta_{i+ i_p}[S^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } |
---|
| 589 | -{\:}_i^k\mathbb{R}_{i_p}^{k_p} |
---|
| 590 | \frac{ \delta_{k+k_p} [S^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } |
---|
| 591 | \right). |
---|
| 592 | \end{equation} |
---|
| 593 | This is symmetrical in $T $ and $S$, so exactly the same term arises |
---|
| 594 | from the discretization of this triad's contribution towards the |
---|
[3218] | 595 | RHS of \eqref{eq:triad:iso_property3}. |
---|
[2282] | 596 | \end{description} |
---|
[3218] | 597 | \subsection{Treatment of the triads at the boundaries}\label{sec:triad:iso_bdry} |
---|
| 598 | The triad slope can only be defined where both the grid boxes centred at |
---|
| 599 | the end of the arms exist. Triads that would poke up |
---|
| 600 | through the upper ocean surface into the atmosphere, or down into the |
---|
| 601 | ocean floor, must be masked out. See Fig.~\ref{fig:triad:bdry_triads}. Surface layer triads |
---|
| 602 | $\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and |
---|
| 603 | $\triad{i+1}{1}{R}{-1/2}{-1/2}$ (blue) that require density to be |
---|
| 604 | specified above the ocean surface are masked (Fig.~\ref{fig:triad:bdry_triads}a): this ensures that lateral |
---|
| 605 | tracer gradients produce no flux through the ocean surface. However, |
---|
| 606 | to prevent surface noise, it is customary to retain the $_{11}$ contributions towards |
---|
| 607 | the lateral triad fluxes $\triad[u]{i}{1}{F}{1/2}{-1/2}$ and |
---|
| 608 | $\triad[u]{i+1}{1}{F}{-1/2}{-1/2}$; this drives diapycnal tracer |
---|
| 609 | fluxes. Similar comments apply to triads that would intersect the |
---|
| 610 | ocean floor (Fig.~\ref{fig:triad:bdry_triads}b). Note that both near bottom |
---|
| 611 | triad slopes $\triad{i}{k}{R}{1/2}{1/2}$ and |
---|
| 612 | $\triad{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the $i,k+1$ |
---|
| 613 | or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point is |
---|
| 614 | masked. The associated lateral fluxes (grey-black dashed line) are |
---|
| 615 | masked if \nlv{ln\_botmix\_grif=.false.}, but left unmasked, |
---|
| 616 | giving bottom mixing, if \nlv{ln\_botmix\_grif=.true.}. |
---|
[2282] | 617 | |
---|
[3218] | 618 | The default option \nlv{ln\_botmix\_grif=.false.} is suitable when the |
---|
| 619 | bbl mixing option is enabled (\key{trabbl}, with \nlv{nn\_bbl\_ldf=1}), |
---|
| 620 | or for simple idealized problems. For setups with topography without |
---|
| 621 | bbl mixing, \nlv{ln\_botmix\_grif=.true.} may be necessary. |
---|
| 622 | % >>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 623 | \begin{figure}[h] \begin{center} |
---|
| 624 | \includegraphics[width=0.60\textwidth]{./TexFiles/Figures/Fig_bdry_triads} |
---|
| 625 | \caption{ \label{fig:triad:bdry_triads} |
---|
| 626 | (a) Uppermost model layer $k=1$ with $i,1$ and $i+1,1$ tracer |
---|
| 627 | points (black dots), and $i+1/2,1$ $u$-point (blue square). Triad |
---|
| 628 | slopes $\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and $\triad{i+1}{1}{R}{-1/2}{-1/2}$ |
---|
| 629 | (blue) poking through the ocean surface are masked (faded in |
---|
| 630 | figure). However, the lateral $_{11}$ contributions towards |
---|
| 631 | $\triad[u]{i}{1}{F}{1/2}{-1/2}$ and $\triad[u]{i+1}{1}{F}{-1/2}{-1/2}$ |
---|
| 632 | (yellow line) are still applied, giving diapycnal diffusive |
---|
| 633 | fluxes.\\ |
---|
| 634 | (b) Both near bottom triad slopes $\triad{i}{k}{R}{1/2}{1/2}$ and |
---|
| 635 | $\triad{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the $i,k+1$ |
---|
| 636 | or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point |
---|
| 637 | is masked. The associated lateral fluxes (grey-black dashed |
---|
| 638 | line) are masked if \smnlv{ln\_botmix\_grif=.false.}, but left |
---|
| 639 | unmasked, giving bottom mixing, if \smnlv{ln\_botmix\_grif=.true.}} |
---|
| 640 | \end{center} \end{figure} |
---|
| 641 | % >>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 642 | \subsection{ Limiting of the slopes within the interior}\label{sec:triad:limit} |
---|
| 643 | As discussed in \S\ref{LDF_slp_iso}, iso-neutral slopes relative to |
---|
| 644 | geopotentials must be bounded everywhere, both for consistency with the small-slope |
---|
| 645 | approximation and for numerical stability \citep{Cox1987, |
---|
| 646 | Griffies_Bk04}. The bound chosen in \NEMO is applied to each |
---|
| 647 | component of the slope separately and has a value of $1/100$ in the ocean interior. |
---|
| 648 | %, ramping linearly down above 70~m depth to zero at the surface |
---|
| 649 | It is of course relevant to the iso-neutral slopes $\tilde{r}_i=r_i+\sigma_i$ relative to |
---|
| 650 | geopotentials (here the $\sigma_i$ are the slopes of the coordinate surfaces relative to |
---|
| 651 | geopotentials) \eqref{Eq_PE_slopes_eiv} rather than the slope $r_i$ relative to coordinate |
---|
| 652 | surfaces, so we require |
---|
| 653 | \begin{equation*} |
---|
| 654 | |\tilde{r}_i|\leq \tilde{r}_\mathrm{max}=0.01. |
---|
| 655 | \end{equation*} |
---|
| 656 | and then recalculate the slopes $r_i$ relative to coordinates. |
---|
| 657 | Each individual triad slope |
---|
| 658 | \begin{equation} |
---|
| 659 | \label{eq:triad:Rtilde} |
---|
| 660 | _i^k\tilde{\mathbb{R}}_{i_p}^{k_p} = {}_i^k\mathbb{R}_{i_p}^{k_p} + \frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}} |
---|
| 661 | \end{equation} |
---|
| 662 | is limited like this and then the corresponding |
---|
| 663 | $_i^k\mathbb{R}_{i_p}^{k_p} $ are recalculated and combined to form the fluxes. |
---|
| 664 | Note that where the slopes have been limited, there is now a non-zero |
---|
| 665 | iso-neutral density flux that drives dianeutral mixing. In particular this iso-neutral density flux |
---|
| 666 | is always downwards, and so acts to reduce gravitational potential energy. |
---|
| 667 | \subsection{Tapering within the surface mixed layer} |
---|
| 668 | Additional tapering of the iso-neutral fluxes is necessary within the |
---|
| 669 | surface mixed layer. When the Griffies triads are used, we offer two |
---|
| 670 | options for this. |
---|
| 671 | \subsubsection{Linear slope tapering within the surface mixed layer}\label{sec:triad:lintaper} |
---|
| 672 | This is the option activated by the default choice |
---|
| 673 | \nlv{ln\_triad\_iso=.false.}. Slopes $\tilde{r}_i$ relative to |
---|
| 674 | geopotentials are tapered linearly from their value immediately below the mixed layer to zero at the |
---|
| 675 | surface, as described in option (c) of Fig.~\ref{Fig_eiv_slp}, to values |
---|
| 676 | \begin{subequations} |
---|
| 677 | \begin{equation} |
---|
| 678 | \label{eq:triad:rmtilde} |
---|
| 679 | \rMLt = |
---|
| 680 | -\frac{z}{h}\left.\tilde{r}_i\right|_{z=-h}\quad \text{ for } z>-h, |
---|
| 681 | \end{equation} |
---|
| 682 | and then the $r_i$ relative to vertical coordinate surfaces are appropriately |
---|
| 683 | adjusted to |
---|
| 684 | \begin{equation} |
---|
| 685 | \label{eq:triad:rm} |
---|
| 686 | \rML =\rMLt -\sigma_i \quad \text{ for } z>-h. |
---|
| 687 | \end{equation} |
---|
| 688 | \end{subequations} |
---|
| 689 | Thus the diffusion operator within the mixed layer is given by: |
---|
| 690 | \begin{equation} \label{eq:triad:iso_tensor_ML} |
---|
| 691 | D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad |
---|
| 692 | \mbox{with}\quad \;\;\Re =\left( {{\begin{array}{*{20}c} |
---|
| 693 | 1 \hfill & 0 \hfill & {-\rML[1]}\hfill \\ |
---|
| 694 | 0 \hfill & 1 \hfill & {-\rML[2]} \hfill \\ |
---|
| 695 | {-\rML[1]}\hfill & {-\rML[2]} \hfill & {\rML[1]^2+\rML[2]^2} \hfill |
---|
| 696 | \end{array} }} \right) |
---|
| 697 | \end{equation} |
---|
[2282] | 698 | |
---|
[3218] | 699 | This slope tapering gives a natural connection between tracer in the |
---|
| 700 | mixed-layer and in isopycnal layers immediately below, in the |
---|
| 701 | thermocline. It is consistent with the way the $\tilde{r}_i$ are |
---|
| 702 | tapered within the mixed layer (see \S\ref{sec:triad:taperskew} below) |
---|
| 703 | so as to ensure a uniform GM eddy-induced velocity throughout the |
---|
| 704 | mixed layer. However, it gives a downwards density flux and so acts so |
---|
| 705 | as to reduce potential energy in the same way as does the slope |
---|
| 706 | limiting discussed above in \S\ref{sec:triad:limit}. |
---|
| 707 | |
---|
| 708 | As in \S\ref{sec:triad:limit} above, the tapering |
---|
| 709 | \eqref{eq:triad:rmtilde} is applied separately to each triad |
---|
| 710 | $_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}$, and the |
---|
| 711 | $_i^k\mathbb{R}_{i_p}^{k_p}$ adjusted. For clarity, we assume |
---|
| 712 | $z$-coordinates in the following; the conversion from |
---|
| 713 | $\mathbb{R}$ to $\tilde{\mathbb{R}}$ and back to $\mathbb{R}$ follows exactly as described |
---|
| 714 | above by \eqref{eq:triad:Rtilde}. |
---|
| 715 | \begin{enumerate} |
---|
| 716 | \item Mixed-layer depth is defined so as to avoid including regions of weak |
---|
| 717 | vertical stratification in the slope definition. |
---|
| 718 | At each $i,j$ (simplified to $i$ in |
---|
| 719 | Fig.~\ref{fig:triad:MLB_triad}), we define the mixed-layer by setting |
---|
| 720 | the vertical index of the tracer point immediately below the mixed |
---|
| 721 | layer, $k_\mathrm{ML}$, as the maximum $k$ (shallowest tracer point) |
---|
| 722 | such that the potential density |
---|
| 723 | ${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$, where $i,k_{10}$ is |
---|
| 724 | the tracer gridbox within which the depth reaches 10~m. See the left |
---|
| 725 | side of Fig.~\ref{fig:triad:MLB_triad}. We use the $k_{10}$-gridbox |
---|
| 726 | instead of the surface gridbox to avoid problems e.g.\ with thin |
---|
| 727 | daytime mixed-layers. Currently we use the same |
---|
| 728 | $\Delta\rho_c=0.01\;\mathrm{kg\:m^{-3}}$ for ML triad tapering as is |
---|
| 729 | used to output the diagnosed mixed-layer depth |
---|
| 730 | $h_\mathrm{ML}=|z_{W}|_{k_\mathrm{ML}+1/2}$, the depth of the $w$-point |
---|
| 731 | above the $i,k_\mathrm{ML}$ tracer point. |
---|
| 732 | |
---|
| 733 | \item We define `basal' triad slopes |
---|
| 734 | ${\:}_i{\mathbb{R}_\mathrm{base}}_{\,i_p}^{k_p}$ as the slopes |
---|
| 735 | of those triads whose vertical `arms' go down from the |
---|
| 736 | $i,k_\mathrm{ML}$ tracer point to the $i,k_\mathrm{ML}-1$ tracer point |
---|
| 737 | below. This is to ensure that the vertical density gradients |
---|
| 738 | associated with these basal triad slopes |
---|
| 739 | ${\:}_i{\mathbb{R}_\mathrm{base}}_{\,i_p}^{k_p}$ are |
---|
| 740 | representative of the thermocline. The four basal triads defined in the bottom part |
---|
| 741 | of Fig.~\ref{fig:triad:MLB_triad} are then |
---|
| 742 | \begin{align} |
---|
| 743 | {\:}_i{\mathbb{R}_\mathrm{base}}_{\,i_p}^{k_p} &= |
---|
| 744 | {\:}^{k_\mathrm{ML}-k_p-1/2}_i{\mathbb{R}_\mathrm{base}}_{\,i_p}^{k_p}, \label{eq:triad:Rbase} |
---|
| 745 | \\ |
---|
| 746 | \intertext{with e.g.\ the green triad} |
---|
| 747 | {\:}_i{\mathbb{R}_\mathrm{base}}_{1/2}^{-1/2}&= |
---|
| 748 | {\:}^{k_\mathrm{ML}}_i{\mathbb{R}_\mathrm{base}}_{\,1/2}^{-1/2}. \notag |
---|
| 749 | \end{align} |
---|
| 750 | The vertical flux associated with each of these triads passes through the $w$-point |
---|
| 751 | $i,k_\mathrm{ML}-1/2$ lying \emph{below} the $i,k_\mathrm{ML}$ tracer point, |
---|
| 752 | so it is this depth |
---|
| 753 | \begin{equation} |
---|
| 754 | \label{eq:triad:zbase} |
---|
| 755 | {z_\mathrm{base}}_{\,i}={z_{w}}_{k_\mathrm{ML}-1/2} |
---|
| 756 | \end{equation} |
---|
| 757 | (one gridbox deeper than the |
---|
| 758 | diagnosed ML depth $z_\mathrm{ML})$ that sets the $h$ used to taper |
---|
| 759 | the slopes in \eqref{eq:triad:rmtilde}. |
---|
| 760 | \item Finally, we calculate the adjusted triads |
---|
| 761 | ${\:}_i^k{\mathbb{R}_\mathrm{ML}}_{\,i_p}^{k_p}$ within the mixed |
---|
| 762 | layer, by multiplying the appropriate |
---|
| 763 | ${\:}_i{\mathbb{R}_\mathrm{base}}_{\,i_p}^{k_p}$ by the ratio of |
---|
| 764 | the depth of the $w$-point ${z_w}_{k+k_p}$ to ${z_\mathrm{base}}_{\,i}$. For |
---|
| 765 | instance the green triad centred on $i,k$ |
---|
| 766 | \begin{align} |
---|
| 767 | {\:}_i^k{\mathbb{R}_\mathrm{ML}}_{\,1/2}^{-1/2} &= |
---|
| 768 | \frac{{z_w}_{k-1/2}}{{z_\mathrm{base}}_{\,i}}{\:}_i{\mathbb{R}_\mathrm{base}}_{\,1/2}^{-1/2} |
---|
| 769 | \notag \\ |
---|
| 770 | \intertext{and more generally} |
---|
| 771 | {\:}_i^k{\mathbb{R}_\mathrm{ML}}_{\,i_p}^{k_p} &= |
---|
| 772 | \frac{{z_w}_{k+k_p}}{{z_\mathrm{base}}_{\,i}}{\:}_i{\mathbb{R}_\mathrm{base}}_{\,i_p}^{k_p}.\label{eq:triad:RML} |
---|
| 773 | \end{align} |
---|
| 774 | \end{enumerate} |
---|
| 775 | |
---|
| 776 | % >>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 777 | \begin{figure}[h] |
---|
| 778 | \fcapside {\caption{\label{fig:triad:MLB_triad} Definition of |
---|
| 779 | mixed-layer depth and calculation of linearly tapered |
---|
| 780 | triads. The figure shows a water column at a given $i,j$ |
---|
| 781 | (simplified to $i$), with the ocean surface at the top. Tracer points are denoted by |
---|
| 782 | bullets, and black lines the edges of the tracer cells; $k$ |
---|
| 783 | increases upwards. \\ |
---|
| 784 | \hspace{5 em}We define the mixed-layer by setting the vertical index |
---|
| 785 | of the tracer point immediately below the mixed layer, |
---|
| 786 | $k_\mathrm{ML}$, as the maximum $k$ (shallowest tracer point) |
---|
| 787 | such that ${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$, |
---|
| 788 | where $i,k_{10}$ is the tracer gridbox within which the depth |
---|
| 789 | reaches 10~m. We calculate the triad slopes within the mixed |
---|
| 790 | layer by linearly tapering them from zero (at the surface) to |
---|
| 791 | the `basal' slopes, the slopes of the four triads passing through the |
---|
| 792 | $w$-point $i,k_\mathrm{ML}-1/2$ (blue square), |
---|
| 793 | ${\:}_i{\mathbb{R}_\mathrm{base}}_{\,i_p}^{k_p}$. Triads with |
---|
| 794 | different $i_p,k_p$, denoted by different colours, (e.g. the green |
---|
| 795 | triad $i_p=1/2,k_p=-1/2$) are tapered to the appropriate basal triad.}} |
---|
| 796 | {\includegraphics[width=0.60\textwidth]{./TexFiles/Figures/Fig_triad_MLB}} |
---|
| 797 | \end{figure} |
---|
| 798 | % >>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 799 | |
---|
| 800 | \subsubsection{Additional truncation of skew iso-neutral flux components} |
---|
| 801 | The alternative option is activated by setting \nlv{ln\_triad\_iso = |
---|
| 802 | .true.}. This retains the same tapered slope $\rML$ described above for the |
---|
| 803 | calculation of the $_{33}$ term of the iso-neutral diffusion tensor (the |
---|
| 804 | vertical tracer flux driven by vertical tracer gradients), but |
---|
| 805 | replaces the $\rML$ in the skew term by |
---|
| 806 | \begin{equation} |
---|
| 807 | \label{eq:triad:rm*} |
---|
| 808 | \rML^*=\left.\rMLt^2\right/\tilde{r}_i-\sigma_i, |
---|
| 809 | \end{equation} |
---|
| 810 | giving a ML diffusive operator |
---|
| 811 | \begin{equation} \label{eq:triad:iso_tensor_ML2} |
---|
| 812 | D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad |
---|
| 813 | \mbox{with}\quad \;\;\Re =\left( {{\begin{array}{*{20}c} |
---|
| 814 | 1 \hfill & 0 \hfill & {-\rML[1]^*}\hfill \\ |
---|
| 815 | 0 \hfill & 1 \hfill & {-\rML[2]^*} \hfill \\ |
---|
| 816 | {-\rML[1]^*}\hfill & {-\rML[2]^*} \hfill & {\rML[1]^2+\rML[2]^2} \hfill \\ |
---|
| 817 | \end{array} }} \right). |
---|
| 818 | \end{equation} |
---|
| 819 | This operator |
---|
| 820 | \footnote{To ensure good behaviour where horizontal density |
---|
| 821 | gradients are weak, we in fact follow \citet{Gerdes1991} and set |
---|
| 822 | $\rML^*=\mathrm{sgn}(\tilde{r}_i)\min(|\rMLt^2/\tilde{r}_i|,|\tilde{r}_i|)-\sigma_i$.} |
---|
| 823 | then has the property it gives no vertical density flux, and so does |
---|
| 824 | not change the potential energy. |
---|
| 825 | This approach is similar to multiplying the iso-neutral diffusion |
---|
| 826 | coefficient by $\tilde{r}_\mathrm{max}^{-2}\tilde{r}_i^{-2}$ for steep |
---|
| 827 | slopes, as suggested by \citet{Gerdes1991} (see also \citet{Griffies_Bk04}). |
---|
| 828 | Again it is applied separately to each triad $_i^k\mathbb{R}_{i_p}^{k_p}$ |
---|
| 829 | |
---|
| 830 | In practice, this approach gives weak vertical tracer fluxes through |
---|
| 831 | the mixed-layer, as well as vanishing density fluxes. While it is |
---|
| 832 | theoretically advantageous that it does not change the potential |
---|
| 833 | energy, it may give a discontinuity between the |
---|
| 834 | fluxes within the mixed-layer (purely horizontal) and just below (along |
---|
| 835 | iso-neutral surfaces). |
---|
| 836 | % This may give strange looking results, |
---|
| 837 | % particularly where the mixed-layer depth varies strongly laterally. |
---|
[2282] | 838 | % ================================================================ |
---|
| 839 | % Skew flux formulation for Eddy Induced Velocity : |
---|
| 840 | % ================================================================ |
---|
[3218] | 841 | \section{Eddy induced advection and its formulation as a skew flux} |
---|
[2282] | 842 | |
---|
[3218] | 843 | \subsection{The continuous skew flux formulation}\label{sec:triad:continuous-skew-flux} |
---|
| 844 | |
---|
| 845 | When Gent and McWilliams's [1990] diffusion is used (\key{traldf\_eiv} defined), |
---|
| 846 | an additional advection term is added. The associated velocity is the so called |
---|
[2282] | 847 | eddy induced velocity, the formulation of which depends on the slopes of iso- |
---|
[3218] | 848 | neutral surfaces. Contrary to the case of iso-neutral mixing, the slopes used |
---|
| 849 | here are referenced to the geopotential surfaces, $i.e.$ \eqref{Eq_ldfslp_geo} |
---|
[2282] | 850 | is used in $z$-coordinate, and the sum \eqref{Eq_ldfslp_geo} |
---|
[3218] | 851 | + \eqref{Eq_ldfslp_iso} in $z^*$ or $s$-coordinates. |
---|
[2282] | 852 | |
---|
[3218] | 853 | The eddy induced velocity is given by: |
---|
| 854 | \begin{equation} \label{eq:triad:eiv_v} |
---|
[2282] | 855 | \begin{split} |
---|
[3218] | 856 | u^* & = - \frac{1}{e_{3}}\; \partial_k \left( A_{e} \; \tilde{r}_1 \right) \\ |
---|
| 857 | v^* & = - \frac{1}{e_{3}}\; \partial_k \left( A_{e} \; \tilde{r}_2 \right) \\ |
---|
| 858 | w^* & = \frac{1}{e_{1}e_{2}}\; \left\{ \partial_i \left( e_{2} \, A_{e} \; \tilde{r}_1 \right) |
---|
| 859 | + \partial_j \left( e_{1} \, A_{e} \;\tilde{r}_2 \right) \right\} |
---|
[2282] | 860 | \end{split} |
---|
| 861 | \end{equation} |
---|
[3218] | 862 | where $A_{e}$ is the eddy induced velocity coefficient, and $\tilde{r}_1$ and $\tilde{r}_2$ the slopes between the iso-neutral and the geopotential surfaces. |
---|
[2282] | 863 | |
---|
[3218] | 864 | The traditional way to implement this additional advection is to add it to the Eulerian |
---|
| 865 | velocity prior to computing the tracer advection. This allows us to take advantage of |
---|
| 866 | all the advection schemes offered for the tracers (see \S\ref{TRA_adv}) and not just |
---|
| 867 | a $2^{nd}$ order advection scheme. This is particularly useful for passive tracers |
---|
| 868 | where \emph{positivity} of the advection scheme is of paramount importance. |
---|
[2282] | 869 | |
---|
[3218] | 870 | \citet{Griffies_JPO98} introduces another way to implement the eddy induced advection, |
---|
| 871 | the so-called skew form. It is based on a transformation of the advective fluxes |
---|
| 872 | using the non-divergent nature of the eddy induced velocity. |
---|
| 873 | For example in the (\textbf{i},\textbf{k}) plane, the tracer advective |
---|
| 874 | fluxes per unit area in $ijk$ space can be |
---|
[2282] | 875 | transformed as follows: |
---|
| 876 | \begin{flalign*} |
---|
| 877 | \begin{split} |
---|
[3218] | 878 | \textbf{F}_\mathrm{eiv}^T = |
---|
| 879 | \begin{pmatrix} |
---|
[2282] | 880 | {e_{2}\,e_{3}\; u^*} \\ |
---|
| 881 | {e_{1}\,e_{2}\; w^*} \\ |
---|
| 882 | \end{pmatrix} \; T |
---|
| 883 | &= |
---|
[3218] | 884 | \begin{pmatrix} |
---|
| 885 | { - \partial_k \left( e_{2} \, A_{e} \; \tilde{r}_1 \right) \; T \;} \\ |
---|
| 886 | {+ \partial_i \left( e_{2} \, A_{e} \; \tilde{r}_1 \right) \; T \;} \\ |
---|
[2282] | 887 | \end{pmatrix} \\ |
---|
[3218] | 888 | &= |
---|
| 889 | \begin{pmatrix} |
---|
| 890 | { - \partial_k \left( e_{2} \, A_{e} \; \tilde{r}_1 \; T \right) \;} \\ |
---|
| 891 | {+ \partial_i \left( e_{2} \, A_{e} \; \tilde{r}_1 \; T \right) \;} \\ |
---|
| 892 | \end{pmatrix} |
---|
| 893 | + |
---|
| 894 | \begin{pmatrix} |
---|
| 895 | {+ e_{2} \, A_{e} \; \tilde{r}_1 \; \partial_k T} \\ |
---|
| 896 | { - e_{2} \, A_{e} \; \tilde{r}_1 \; \partial_i T} \\ |
---|
| 897 | \end{pmatrix} |
---|
[2282] | 898 | \end{split} |
---|
| 899 | \end{flalign*} |
---|
[3218] | 900 | and since the eddy induced velocity field is non-divergent, we end up with the skew |
---|
| 901 | form of the eddy induced advective fluxes per unit area in $ijk$ space: |
---|
| 902 | \begin{equation} \label{eq:triad:eiv_skew_ijk} |
---|
| 903 | \textbf{F}_\mathrm{eiv}^T = \begin{pmatrix} |
---|
| 904 | {+ e_{2} \, A_{e} \; \tilde{r}_1 \; \partial_k T} \\ |
---|
| 905 | { - e_{2} \, A_{e} \; \tilde{r}_1 \; \partial_i T} \\ |
---|
[2282] | 906 | \end{pmatrix} |
---|
| 907 | \end{equation} |
---|
[3218] | 908 | The total fluxes per unit physical area are then |
---|
| 909 | \begin{equation}\label{eq:triad:eiv_skew_physical} |
---|
| 910 | \begin{split} |
---|
| 911 | f^*_1 & = \frac{1}{e_{3}}\; A_{e} \; \tilde{r}_1 \partial_k T \\ |
---|
| 912 | f^*_2 & = \frac{1}{e_{3}}\; A_{e} \; \tilde{r}_2 \partial_k T \\ |
---|
| 913 | f^*_3 & = -\frac{1}{e_{1}e_{2}}\; A_{e} \left\{ e_{2} \tilde{r}_1 \partial_i T |
---|
| 914 | + e_{1} \tilde{r}_2 \partial_j T \right\}. \\ |
---|
| 915 | \end{split} |
---|
| 916 | \end{equation} |
---|
| 917 | Note that Eq.~ \eqref{eq:triad:eiv_skew_physical} takes the same form whatever the |
---|
| 918 | vertical coordinate, though of course the slopes |
---|
| 919 | $\tilde{r}_i$ are relative to geopotentials. |
---|
| 920 | The tendency associated with eddy induced velocity is then simply the convergence |
---|
| 921 | of the fluxes (\ref{eq:triad:eiv_skew_ijk}, \ref{eq:triad:eiv_skew_physical}), so |
---|
| 922 | \begin{equation} \label{eq:triad:skew_eiv_conv} |
---|
| 923 | \frac{\partial T}{\partial t}= -\frac{1}{e_1 \, e_2 \, e_3 } \left[ |
---|
| 924 | \frac{\partial}{\partial i} \left( e_2 A_{e} \; \tilde{r}_1 \partial_k T\right) |
---|
| 925 | + \frac{\partial}{\partial j} \left( e_1 A_{e} \; |
---|
| 926 | \tilde{r}_2 \partial_k T\right) |
---|
| 927 | - \frac{\partial}{\partial k} A_{e} \left( e_{2} \tilde{r}_1 \partial_i T |
---|
| 928 | + e_{1} \tilde{r}_2 \partial_j T \right) \right] |
---|
| 929 | \end{equation} |
---|
| 930 | It naturally conserves the tracer content, as it is expressed in flux |
---|
| 931 | form. Since it has the same divergence as the advective form it also |
---|
| 932 | preserves the tracer variance. |
---|
[2282] | 933 | |
---|
[3218] | 934 | \subsection{The discrete skew flux formulation} |
---|
| 935 | The skew fluxes in (\ref{eq:triad:eiv_skew_physical}, \ref{eq:triad:eiv_skew_ijk}), like the off-diagonal terms |
---|
| 936 | (\ref{eq:triad:i13c}, \ref{eq:triad:i31c}) of the small angle diffusion tensor, are best |
---|
| 937 | expressed in terms of the triad slopes, as in Fig.~\ref{fig:triad:ISO_triad} |
---|
| 938 | and Eqs~(\ref{eq:triad:i13}, \ref{eq:triad:i31}); but now in terms of the triad slopes |
---|
| 939 | $\tilde{\mathbb{R}}$ relative to geopotentials instead of the |
---|
| 940 | $\mathbb{R}$ relative to coordinate surfaces. The discrete form of |
---|
| 941 | \eqref{eq:triad:eiv_skew_ijk} using the slopes \eqref{eq:triad:R} and |
---|
| 942 | defining $A_e$ at $T$-points is then given by: |
---|
[2282] | 943 | |
---|
| 944 | |
---|
[3218] | 945 | \begin{subequations}\label{eq:triad:allskewflux} |
---|
| 946 | \begin{flalign}\label{eq:triad:vect_skew_flux} |
---|
| 947 | \vect{F}_\mathrm{eiv}(T) &\equiv |
---|
| 948 | \sum_{\substack{i_p,\,k_p}} |
---|
| 949 | \begin{pmatrix} |
---|
| 950 | {_{i+1/2-i_p}^k {\mathbb{S}_u}_{i_p}^{k_p} } (T) \\ |
---|
| 951 | \\ |
---|
| 952 | {_i^{k+1/2-k_p} {\mathbb{S}_w}_{i_p}^{k_p} } (T) \\ |
---|
| 953 | \end{pmatrix}, |
---|
| 954 | \end{flalign} |
---|
| 955 | where the skew flux in the $i$-direction associated with a given |
---|
| 956 | triad is (\ref{eq:triad:latflux-triad}, \ref{eq:triad:triadfluxu}): |
---|
| 957 | \begin{align} |
---|
| 958 | \label{eq:triad:skewfluxu} |
---|
| 959 | _i^k {\mathbb{S}_u}_{i_p}^{k_p} (T) &= + \quarter {A_e}_i^k{ |
---|
| 960 | \:}\frac{{b_u}_{i+i_p}^k}{{e_{1u}}_{\,i + i_p}^{\,k}} |
---|
| 961 | \ {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}} \ |
---|
| 962 | \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }, |
---|
| 963 | \\ |
---|
| 964 | \intertext{and \eqref{eq:triad:triadfluxw} in the $k$-direction, changing the sign |
---|
| 965 | to be consistent with \eqref{eq:triad:eiv_skew_ijk}:} |
---|
| 966 | _i^k {\mathbb{S}_w}_{i_p}^{k_p} (T) |
---|
| 967 | &= -\quarter {A_e}_i^k{\: }\frac{{b_u}_{i+i_p}^k}{{e_{3w}}_{\,i}^{\,k+k_p}} |
---|
| 968 | {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }.\label{eq:triad:skewfluxw} |
---|
| 969 | \end{align} |
---|
| 970 | \end{subequations} |
---|
[2282] | 971 | |
---|
[3218] | 972 | Such a discretisation is consistent with the iso-neutral |
---|
| 973 | operator as it uses the same definition for the slopes. It also |
---|
| 974 | ensures the following two key properties. |
---|
| 975 | \subsubsection{No change in tracer variance} |
---|
| 976 | The discretization conserves tracer variance, $i.e.$ it does not |
---|
| 977 | include a diffusive component but is a `pure' advection term. This can |
---|
| 978 | be seen either from Appendix \ref{Apdx_eiv_skew} or by considering the |
---|
| 979 | fluxes associated with a given triad slope |
---|
| 980 | $_i^k{\mathbb{R}}_{i_p}^{k_p} (T)$. For, following |
---|
| 981 | \S\ref{sec:triad:variance} and \eqref{eq:triad:dvar_iso_i}, the |
---|
| 982 | associated horizontal skew-flux $_i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)$ |
---|
| 983 | drives a net rate of change of variance, summed over the two |
---|
| 984 | $T$-points $i+i_p-\half,k$ and $i+i_p+\half,k$, of |
---|
| 985 | \begin{equation} |
---|
| 986 | \label{eq:triad:dvar_eiv_i} |
---|
| 987 | _i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], |
---|
| 988 | \end{equation} |
---|
| 989 | while the associated vertical skew-flux gives a variance change summed over the |
---|
| 990 | $T$-points $i,k+k_p-\half$ (above) and $i,k+k_p+\half$ (below) of |
---|
| 991 | \begin{equation} |
---|
| 992 | \label{eq:triad:dvar_eiv_k} |
---|
| 993 | _i^k{\mathbb{S}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i]. |
---|
| 994 | \end{equation} |
---|
| 995 | Inspection of the definitions (\ref{eq:triad:skewfluxu}, \ref{eq:triad:skewfluxw}) |
---|
| 996 | shows that these two variance changes (\ref{eq:triad:dvar_eiv_i}, \ref{eq:triad:dvar_eiv_k}) |
---|
| 997 | sum to zero. Hence the two fluxes associated with each triad make no |
---|
| 998 | net contribution to the variance budget. |
---|
[2282] | 999 | |
---|
[3218] | 1000 | \subsubsection{Reduction in gravitational PE} |
---|
| 1001 | The vertical density flux associated with the vertical skew-flux |
---|
| 1002 | always has the same sign as the vertical density gradient; thus, so |
---|
| 1003 | long as the fluid is stable (the vertical density gradient is |
---|
| 1004 | negative) the vertical density flux is negative (downward) and hence |
---|
| 1005 | reduces the gravitational PE. |
---|
[2282] | 1006 | |
---|
[3218] | 1007 | For the change in gravitational PE driven by the $k$-flux is |
---|
| 1008 | \begin{align} |
---|
| 1009 | \label{eq:triad:vert_densityPE} |
---|
| 1010 | g {e_{3w}}_{\,i}^{\,k+k_p}{\mathbb{S}_w}_{i_p}^{k_p} (\rho) |
---|
| 1011 | &=g {e_{3w}}_{\,i}^{\,k+k_p}\left[-\alpha _i^k {\:}_i^k |
---|
| 1012 | {\mathbb{S}_w}_{i_p}^{k_p} (T) + \beta_i^k {\:}_i^k |
---|
| 1013 | {\mathbb{S}_w}_{i_p}^{k_p} (S) \right]. \notag \\ |
---|
| 1014 | \intertext{Substituting ${\:}_i^k {\mathbb{S}_w}_{i_p}^{k_p}$ from |
---|
| 1015 | \eqref{eq:triad:skewfluxw}, gives} |
---|
| 1016 | % and separating out |
---|
| 1017 | % $\rtriadt{R}=\rtriad{R} + \delta_{i+i_p}[z_T^k]$, |
---|
| 1018 | % gives two terms. The |
---|
| 1019 | % first $\rtriad{R}$ term (the only term for $z$-coordinates) is: |
---|
| 1020 | &=-\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}} |
---|
| 1021 | \frac{ -\alpha _i^k\delta_{i+ i_p}[T^k]+ \beta_i^k\delta_{i+ i_p}[S^k]} { {e_{1u}}_{\,i + i_p}^{\,k} } \notag \\ |
---|
| 1022 | &=+\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k |
---|
| 1023 | \left({_i^k\mathbb{R}_{i_p}^{k_p}}+\frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}\right) {_i^k\mathbb{R}_{i_p}^{k_p}} |
---|
| 1024 | \frac{-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]} {{e_{3w}}_{\,i}^{\,k+k_p}}, |
---|
| 1025 | \end{align} |
---|
| 1026 | using the definition of the triad slope $\rtriad{R}$, |
---|
| 1027 | \eqref{eq:triad:R} to express $-\alpha _i^k\delta_{i+ i_p}[T^k]+ |
---|
| 1028 | \beta_i^k\delta_{i+ i_p}[S^k]$ in terms of $-\alpha_i^k \delta_{k+ |
---|
| 1029 | k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]$. |
---|
| 1030 | |
---|
| 1031 | Where the coordinates slope, the $i$-flux gives a PE change |
---|
| 1032 | \begin{multline} |
---|
| 1033 | \label{eq:triad:lat_densityPE} |
---|
| 1034 | g \delta_{i+i_p}[z_T^k] |
---|
| 1035 | \left[ |
---|
| 1036 | -\alpha _i^k {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (T) + \beta_i^k {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (S) |
---|
| 1037 | \right] \\ |
---|
| 1038 | = +\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k |
---|
| 1039 | \frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}} |
---|
| 1040 | \left({_i^k\mathbb{R}_{i_p}^{k_p}}+\frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}\right) |
---|
| 1041 | \frac{-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]} {{e_{3w}}_{\,i}^{\,k+k_p}}, |
---|
| 1042 | \end{multline} |
---|
| 1043 | (using \eqref{eq:triad:skewfluxu}) and so the total PE change |
---|
| 1044 | \eqref{eq:triad:vert_densityPE} + \eqref{eq:triad:lat_densityPE} associated with the triad fluxes is |
---|
| 1045 | \begin{multline} |
---|
| 1046 | \label{eq:triad:tot_densityPE} |
---|
| 1047 | g{e_{3w}}_{\,i}^{\,k+k_p}{\mathbb{S}_w}_{i_p}^{k_p} (\rho) + |
---|
| 1048 | g\delta_{i+i_p}[z_T^k] {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (\rho) \\ |
---|
| 1049 | = +\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k |
---|
| 1050 | \left({_i^k\mathbb{R}_{i_p}^{k_p}}+\frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}\right)^2 |
---|
| 1051 | \frac{-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]} {{e_{3w}}_{\,i}^{\,k+k_p}}. |
---|
| 1052 | \end{multline} |
---|
| 1053 | Where the fluid is stable, with $-\alpha_i^k \delta_{k+ k_p}[T^i]+ |
---|
| 1054 | \beta_i^k\delta_{k+ k_p}[S^i]<0$, this PE change is negative. |
---|
| 1055 | |
---|
| 1056 | \subsection{Treatment of the triads at the boundaries}\label{sec:triad:skew_bdry} |
---|
| 1057 | Triad slopes \rtriadt{R} used for the calculation of the eddy-induced skew-fluxes |
---|
| 1058 | are masked at the boundaries in exactly the same way as are the triad |
---|
| 1059 | slopes \rtriad{R} used for the iso-neutral diffusive fluxes, as |
---|
| 1060 | described in \S\ref{sec:triad:iso_bdry} and |
---|
| 1061 | Fig.~\ref{fig:triad:bdry_triads}. Thus surface layer triads |
---|
| 1062 | $\triadt{i}{1}{R}{1/2}{-1/2}$ and $\triadt{i+1}{1}{R}{-1/2}{-1/2}$ are |
---|
| 1063 | masked, and both near bottom triad slopes $\triadt{i}{k}{R}{1/2}{1/2}$ |
---|
| 1064 | and $\triadt{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the |
---|
| 1065 | $i,k+1$ or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ |
---|
| 1066 | $u$-point is masked. The namelist parameter \nlv{ln\_botmix\_grif} has |
---|
| 1067 | no effect on the eddy-induced skew-fluxes. |
---|
| 1068 | |
---|
| 1069 | \subsection{ Limiting of the slopes within the interior}\label{sec:triad:limitskew} |
---|
| 1070 | Presently, the iso-neutral slopes $\tilde{r}_i$ relative |
---|
| 1071 | to geopotentials are limited to be less than $1/100$, exactly as in |
---|
| 1072 | calculating the iso-neutral diffusion, \S \ref{sec:triad:limit}. Each |
---|
| 1073 | individual triad \rtriadt{R} is so limited. |
---|
| 1074 | |
---|
| 1075 | \subsection{Tapering within the surface mixed layer}\label{sec:triad:taperskew} |
---|
| 1076 | The slopes $\tilde{r}_i$ relative to |
---|
| 1077 | geopotentials (and thus the individual triads \rtriadt{R}) are always tapered linearly from their value immediately below the mixed layer to zero at the |
---|
| 1078 | surface \eqref{eq:triad:rmtilde}, as described in \S\ref{sec:triad:lintaper}. This is |
---|
| 1079 | option (c) of Fig.~\ref{Fig_eiv_slp}. This linear tapering for the |
---|
| 1080 | slopes used to calculate the eddy-induced fluxes is |
---|
| 1081 | unaffected by the value of \nlv{ln\_triad\_iso}. |
---|
| 1082 | |
---|
| 1083 | The justification for this linear slope tapering is that, for $A_e$ |
---|
| 1084 | that is constant or varies only in the horizontal (the most commonly |
---|
| 1085 | used options in \NEMO: see \S\ref{LDF_coef}), it is |
---|
| 1086 | equivalent to a horizontal eiv (eddy-induced velocity) that is uniform |
---|
| 1087 | within the mixed layer \eqref{eq:triad:eiv_v}. This ensures that the |
---|
| 1088 | eiv velocities do not restratify the mixed layer \citep{Treguier1997, |
---|
| 1089 | Danabasoglu_al_2008}. Equivantly, in terms |
---|
| 1090 | of the skew-flux formulation we use here, the |
---|
| 1091 | linear slope tapering within the mixed-layer gives a linearly varying |
---|
| 1092 | vertical flux, and so a tracer convergence uniform in depth (the |
---|
| 1093 | horizontal flux convergence is relatively insignificant within the mixed-layer). |
---|
| 1094 | |
---|
| 1095 | \subsection{Streamfunction diagnostics}\label{sec:triad:sfdiag} |
---|
| 1096 | Where the namelist parameter \nlv{ln\_botmix\_grif=.true.}, diagnosed |
---|
| 1097 | mean eddy-induced velocities are output. Each time step, |
---|
| 1098 | streamfunctions are calculated in the $i$-$k$ and $j$-$k$ planes at |
---|
| 1099 | $uw$ (integer +1/2 $i$, integer $j$, integer +1/2 $k$) and $vw$ |
---|
| 1100 | (integer $i$, integer +1/2 $j$, integer +1/2 $k$) points (see Table |
---|
| 1101 | \ref{Tab_cell}) respectively. We follow \citep{Griffies_Bk04} and |
---|
| 1102 | calculate the streamfunction at a given $uw$-point from the |
---|
| 1103 | surrounding four triads according to: |
---|
| 1104 | \begin{equation} |
---|
| 1105 | \label{eq:triad:sfdiagi} |
---|
| 1106 | {\psi_{[i]}}_{i+1/2}^{k+1/2}={\quarter}\sum_{\substack{i_p,\,k_p}} |
---|
| 1107 | {A_e}_{i+1/2-i_p}^{k+1/2-k_p}\:\triadd{i+1/2-i_p}{k+1/2-k_p}{R}{i_p}{k_p} |
---|
| 1108 | \end{equation} |
---|
| 1109 | |
---|
[2282] | 1110 | \newpage %force an empty line |
---|
| 1111 | % ================================================================ |
---|
| 1112 | % Discrete Invariants of the skew flux formulation |
---|
| 1113 | % ================================================================ |
---|
| 1114 | \subsection{Discrete Invariants of the skew flux formulation} |
---|
| 1115 | \label{Apdx_eiv_skew} |
---|
| 1116 | |
---|
| 1117 | |
---|
[3218] | 1118 | Demonstration for the conservation of the tracer variance in the (\textbf{i},\textbf{j}) plane. |
---|
[2282] | 1119 | |
---|
| 1120 | This have to be moved in an Appendix. |
---|
| 1121 | |
---|
| 1122 | The continuous property to be demonstrated is : |
---|
| 1123 | \begin{align*} |
---|
[3218] | 1124 | \int_D \nabla \cdot \textbf{F}_\mathrm{eiv}(T) \; T \;dv \equiv 0 |
---|
[2282] | 1125 | \end{align*} |
---|
[3218] | 1126 | The discrete form of its left hand side is obtained using \eqref{eq:triad:allskewflux} |
---|
[2282] | 1127 | \begin{align*} |
---|
| 1128 | \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}} \Biggl\{ \;\; |
---|
[3218] | 1129 | \delta_i &\left[ |
---|
| 1130 | {e_{2u}}_{i+i_p+1/2}^{k} \;\ \ {A_{e}}_{i+i_p+1/2}^{k} |
---|
| 1131 | \ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} } \quad \delta_{k+k_p}[T_{i+i_p+1/2}] |
---|
[2282] | 1132 | \right] \; T_i^k \\ |
---|
[3218] | 1133 | - \delta_k &\left[ |
---|
| 1134 | {e_{2u}}_i^{k+k_p+1/2} \ \ {A_{e}}_i^{k+k_p+1/2} |
---|
| 1135 | \ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} } \ \ \delta_{i+i_p}[T^{k+k_p+1/2}] |
---|
| 1136 | \right] \; T_i^k \ \Biggr\} |
---|
[2282] | 1137 | \end{align*} |
---|
| 1138 | apply the adjoint of delta operator, it becomes |
---|
| 1139 | \begin{align*} |
---|
| 1140 | \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}} \Biggl\{ \;\; |
---|
[3218] | 1141 | &\left( |
---|
| 1142 | {e_{2u}}_{i+i_p+1/2}^{k} \;\ \ {A_{e}}_{i+i_p+1/2}^{k} |
---|
| 1143 | \ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} } \quad \delta_{k+k_p}[T_{i+i_p+1/2}] |
---|
[2282] | 1144 | \right) \; \delta_{i+1/2}[T^{k}] \\ |
---|
[3218] | 1145 | - &\left( |
---|
| 1146 | {e_{2u}}_i^{k+k_p+1/2} \ \ {A_{e}}_i^{k+k_p+1/2} |
---|
| 1147 | \ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} } \ \ \delta_{i+i_p}[T^{k+k_p+1/2}] |
---|
| 1148 | \right) \; \delta_{k+1/2}[T_{i}] \ \Biggr\} |
---|
[2282] | 1149 | \end{align*} |
---|
| 1150 | Expending the summation on $i_p$ and $k_p$, it becomes: |
---|
| 1151 | \begin{align*} |
---|
[3218] | 1152 | \begin{matrix} |
---|
| 1153 | &\sum\limits_{i,k} \Bigl\{ |
---|
| 1154 | &+{e_{2u}}_{i+1}^{k} &{A_{e}}_{i+1 }^{k} |
---|
[2282] | 1155 | &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{-1/2}} &\delta_{k-1/2}[T_{i+1}] &\delta_{i+1/2}[T^{k}] &\\ |
---|
[3218] | 1156 | &&+{e_{2u}}_i^{k\ \ \ \:} &{A_{e}}_{i}^{k\ \ \ \:} |
---|
[2282] | 1157 | &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{-1/2}} &\delta_{k-1/2}[T_{i\ \ \ \;}] &\delta_{i+1/2}[T^{k}] &\\ |
---|
[3218] | 1158 | &&+{e_{2u}}_{i+1}^{k} &{A_{e}}_{i+1 }^{k} |
---|
[2282] | 1159 | &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{+1/2}} &\delta_{k+1/2}[T_{i+1}] &\delta_{i+1/2}[T^{k}] &\\ |
---|
[3218] | 1160 | &&+{e_{2u}}_i^{k\ \ \ \:} &{A_{e}}_{i}^{k\ \ \ \:} |
---|
[2282] | 1161 | &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{+1/2}} &\delta_{k+1/2}[T_{i\ \ \ \;}] &\delta_{i+1/2}[T^{k}] &\\ |
---|
| 1162 | % |
---|
[3218] | 1163 | &&-{e_{2u}}_i^{k+1} &{A_{e}}_i^{k+1} |
---|
| 1164 | &{_i^{k+1} \mathbb{R}_{-1/2}^{- 1/2}} &\delta_{i-1/2}[T^{k+1}] &\delta_{k+1/2}[T_{i}] &\\ |
---|
| 1165 | &&-{e_{2u}}_i^{k\ \ \ \:} &{A_{e}}_i^{k\ \ \ \:} |
---|
[2282] | 1166 | &{\ \ \;_i^k \mathbb{R}_{-1/2}^{+1/2}} &\delta_{i-1/2}[T^{k\ \ \ \:}] &\delta_{k+1/2}[T_{i}] &\\ |
---|
[3218] | 1167 | &&-{e_{2u}}_i^{k+1 } &{A_{e}}_i^{k+1} |
---|
| 1168 | &{_i^{k+1} \mathbb{R}_{+1/2}^{- 1/2}} &\delta_{i+1/2}[T^{k+1}] &\delta_{k+1/2}[T_{i}] &\\ |
---|
| 1169 | &&-{e_{2u}}_i^{k\ \ \ \:} &{A_{e}}_i^{k\ \ \ \:} |
---|
| 1170 | &{\ \ \;_i^k \mathbb{R}_{+1/2}^{+1/2}} &\delta_{i+1/2}[T^{k\ \ \ \:}] &\delta_{k+1/2}[T_{i}] |
---|
[2282] | 1171 | &\Bigr\} \\ |
---|
[3218] | 1172 | \end{matrix} |
---|
[2282] | 1173 | \end{align*} |
---|
[3218] | 1174 | The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{+1/2}}$ are the |
---|
| 1175 | same but of opposite signs, they cancel out. |
---|
| 1176 | Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{-1/2}}$. |
---|
| 1177 | The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{-1/2}}$ are the |
---|
| 1178 | same but both of opposite signs and shifted by 1 in $k$ direction. When summing over $k$ |
---|
| 1179 | they cancel out with the neighbouring grid points. |
---|
| 1180 | Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{+1/2}}$ in the |
---|
| 1181 | $i$ direction. Therefore the sum over the domain is zero, $i.e.$ the variance of the |
---|
[2282] | 1182 | tracer is preserved by the discretisation of the skew fluxes. |
---|
| 1183 | |
---|
[3218] | 1184 | %%% Local Variables: |
---|
| 1185 | %%% TeX-master: "../../NEMO_book-luatex.tex" |
---|
| 1186 | %%% End: |
---|