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