1 | \documentclass[NEMO_book]{subfiles} |
---|
2 | \begin{document} |
---|
3 | % ================================================================ |
---|
4 | % Chapter Vertical Ocean Physics (ZDF) |
---|
5 | % ================================================================ |
---|
6 | \chapter{Vertical Ocean Physics (ZDF)} |
---|
7 | \label{ZDF} |
---|
8 | \minitoc |
---|
9 | |
---|
10 | %gm% Add here a small introduction to ZDF and naming of the different physics (similar to what have been written for TRA and DYN. |
---|
11 | |
---|
12 | |
---|
13 | \newpage |
---|
14 | $\ $\newline % force a new ligne |
---|
15 | |
---|
16 | |
---|
17 | % ================================================================ |
---|
18 | % Vertical Mixing |
---|
19 | % ================================================================ |
---|
20 | \section{Vertical Mixing} |
---|
21 | \label{ZDF_zdf} |
---|
22 | |
---|
23 | The discrete form of the ocean subgrid scale physics has been presented in |
---|
24 | \S\ref{TRA_zdf} and \S\ref{DYN_zdf}. At the surface and bottom boundaries, |
---|
25 | the turbulent fluxes of momentum, heat and salt have to be defined. At the |
---|
26 | surface they are prescribed from the surface forcing (see Chap.~\ref{SBC}), |
---|
27 | while at the bottom they are set to zero for heat and salt, unless a geothermal |
---|
28 | flux forcing is prescribed as a bottom boundary condition ($i.e.$ \key{trabbl} |
---|
29 | defined, see \S\ref{TRA_bbc}), and specified through a bottom friction |
---|
30 | parameterisation for momentum (see \S\ref{ZDF_bfr}). |
---|
31 | |
---|
32 | In this section we briefly discuss the various choices offered to compute |
---|
33 | the vertical eddy viscosity and diffusivity coefficients, $A_u^{vm}$ , |
---|
34 | $A_v^{vm}$ and $A^{vT}$ ($A^{vS}$), defined at $uw$-, $vw$- and $w$- |
---|
35 | points, respectively (see \S\ref{TRA_zdf} and \S\ref{DYN_zdf}). These |
---|
36 | coefficients can be assumed to be either constant, or a function of the local |
---|
37 | Richardson number, or computed from a turbulent closure model (either |
---|
38 | TKE or GLS formulation). The computation of these coefficients is initialized |
---|
39 | in the \mdl{zdfini} module and performed in the \mdl{zdfric}, \mdl{zdftke} or |
---|
40 | \mdl{zdfgls} modules. The trends due to the vertical momentum and tracer |
---|
41 | diffusion, including the surface forcing, are computed and added to the |
---|
42 | general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively. |
---|
43 | These trends can be computed using either a forward time stepping scheme |
---|
44 | (namelist parameter \np{ln\_zdfexp}=true) or a backward time stepping |
---|
45 | scheme (\np{ln\_zdfexp}=false) depending on the magnitude of the mixing |
---|
46 | coefficients, and thus of the formulation used (see \S\ref{STP}). |
---|
47 | |
---|
48 | % ------------------------------------------------------------------------------------------------------------- |
---|
49 | % Constant |
---|
50 | % ------------------------------------------------------------------------------------------------------------- |
---|
51 | \subsection{Constant (\protect\key{zdfcst})} |
---|
52 | \label{ZDF_cst} |
---|
53 | %--------------------------------------------namzdf--------------------------------------------------------- |
---|
54 | \namdisplay{namzdf} |
---|
55 | %-------------------------------------------------------------------------------------------------------------- |
---|
56 | |
---|
57 | Options are defined through the \ngn{namzdf} namelist variables. |
---|
58 | When \key{zdfcst} is defined, the momentum and tracer vertical eddy coefficients |
---|
59 | are set to constant values over the whole ocean. This is the crudest way to define |
---|
60 | the vertical ocean physics. It is recommended that this option is only used in |
---|
61 | process studies, not in basin scale simulations. Typical values used in this case are: |
---|
62 | \begin{align*} |
---|
63 | A_u^{vm} = A_v^{vm} &= 1.2\ 10^{-4}~m^2.s^{-1} \\ |
---|
64 | A^{vT} = A^{vS} &= 1.2\ 10^{-5}~m^2.s^{-1} |
---|
65 | \end{align*} |
---|
66 | |
---|
67 | These values are set through the \np{rn\_avm0} and \np{rn\_avt0} namelist parameters. |
---|
68 | In all cases, do not use values smaller that those associated with the molecular |
---|
69 | viscosity and diffusivity, that is $\sim10^{-6}~m^2.s^{-1}$ for momentum, |
---|
70 | $\sim10^{-7}~m^2.s^{-1}$ for temperature and $\sim10^{-9}~m^2.s^{-1}$ for salinity. |
---|
71 | |
---|
72 | |
---|
73 | % ------------------------------------------------------------------------------------------------------------- |
---|
74 | % Richardson Number Dependent |
---|
75 | % ------------------------------------------------------------------------------------------------------------- |
---|
76 | \subsection{Richardson Number Dependent (\protect\key{zdfric})} |
---|
77 | \label{ZDF_ric} |
---|
78 | |
---|
79 | %--------------------------------------------namric--------------------------------------------------------- |
---|
80 | \namdisplay{namzdf_ric} |
---|
81 | %-------------------------------------------------------------------------------------------------------------- |
---|
82 | |
---|
83 | When \key{zdfric} is defined, a local Richardson number dependent formulation |
---|
84 | for the vertical momentum and tracer eddy coefficients is set through the \ngn{namzdf\_ric} |
---|
85 | namelist variables.The vertical mixing |
---|
86 | coefficients are diagnosed from the large scale variables computed by the model. |
---|
87 | \textit{In situ} measurements have been used to link vertical turbulent activity to |
---|
88 | large scale ocean structures. The hypothesis of a mixing mainly maintained by the |
---|
89 | growth of Kelvin-Helmholtz like instabilities leads to a dependency between the |
---|
90 | vertical eddy coefficients and the local Richardson number ($i.e.$ the |
---|
91 | ratio of stratification to vertical shear). Following \citet{Pacanowski_Philander_JPO81}, the following |
---|
92 | formulation has been implemented: |
---|
93 | \begin{equation} \label{Eq_zdfric} |
---|
94 | \left\{ \begin{aligned} |
---|
95 | A^{vT} &= \frac {A_{ric}^{vT}}{\left( 1+a \; Ri \right)^n} + A_b^{vT} \\ |
---|
96 | A^{vm} &= \frac{A^{vT} }{\left( 1+ a \;Ri \right) } + A_b^{vm} |
---|
97 | \end{aligned} \right. |
---|
98 | \end{equation} |
---|
99 | where $Ri = N^2 / \left(\partial_z \textbf{U}_h \right)^2$ is the local Richardson |
---|
100 | number, $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}), |
---|
101 | $A_b^{vT} $ and $A_b^{vm}$ are the constant background values set as in the |
---|
102 | constant case (see \S\ref{ZDF_cst}), and $A_{ric}^{vT} = 10^{-4}~m^2.s^{-1}$ |
---|
103 | is the maximum value that can be reached by the coefficient when $Ri\leq 0$, |
---|
104 | $a=5$ and $n=2$. The last three values can be modified by setting the |
---|
105 | \np{rn\_avmri}, \np{rn\_alp} and \np{nn\_ric} namelist parameters, respectively. |
---|
106 | |
---|
107 | A simple mixing-layer model to transfer and dissipate the atmospheric |
---|
108 | forcings (wind-stress and buoyancy fluxes) can be activated setting |
---|
109 | the \np{ln\_mldw} =.true. in the namelist. |
---|
110 | |
---|
111 | In this case, the local depth of turbulent wind-mixing or "Ekman depth" |
---|
112 | $h_{e}(x,y,t)$ is evaluated and the vertical eddy coefficients prescribed within this layer. |
---|
113 | |
---|
114 | This depth is assumed proportional to the "depth of frictional influence" that is limited by rotation: |
---|
115 | \begin{equation} |
---|
116 | h_{e} = Ek \frac {u^{*}} {f_{0}} \\ |
---|
117 | \end{equation} |
---|
118 | where, $Ek$ is an empirical parameter, $u^{*}$ is the friction velocity and $f_{0}$ is the Coriolis |
---|
119 | parameter. |
---|
120 | |
---|
121 | In this similarity height relationship, the turbulent friction velocity: |
---|
122 | \begin{equation} |
---|
123 | u^{*} = \sqrt \frac {|\tau|} {\rho_o} \\ |
---|
124 | \end{equation} |
---|
125 | |
---|
126 | is computed from the wind stress vector $|\tau|$ and the reference density $ \rho_o$. |
---|
127 | The final $h_{e}$ is further constrained by the adjustable bounds \np{rn\_mldmin} and \np{rn\_mldmax}. |
---|
128 | Once $h_{e}$ is computed, the vertical eddy coefficients within $h_{e}$ are set to |
---|
129 | the empirical values \np{rn\_wtmix} and \np{rn\_wvmix} \citep{Lermusiaux2001}. |
---|
130 | |
---|
131 | % ------------------------------------------------------------------------------------------------------------- |
---|
132 | % TKE Turbulent Closure Scheme |
---|
133 | % ------------------------------------------------------------------------------------------------------------- |
---|
134 | \subsection{TKE Turbulent Closure Scheme (\protect\key{zdftke})} |
---|
135 | \label{ZDF_tke} |
---|
136 | |
---|
137 | %--------------------------------------------namzdf_tke-------------------------------------------------- |
---|
138 | \namdisplay{namzdf_tke} |
---|
139 | %-------------------------------------------------------------------------------------------------------------- |
---|
140 | |
---|
141 | The vertical eddy viscosity and diffusivity coefficients are computed from a TKE |
---|
142 | turbulent closure model based on a prognostic equation for $\bar{e}$, the turbulent |
---|
143 | kinetic energy, and a closure assumption for the turbulent length scales. This |
---|
144 | turbulent closure model has been developed by \citet{Bougeault1989} in the |
---|
145 | atmospheric case, adapted by \citet{Gaspar1990} for the oceanic case, and |
---|
146 | embedded in OPA, the ancestor of NEMO, by \citet{Blanke1993} for equatorial Atlantic |
---|
147 | simulations. Since then, significant modifications have been introduced by |
---|
148 | \citet{Madec1998} in both the implementation and the formulation of the mixing |
---|
149 | length scale. The time evolution of $\bar{e}$ is the result of the production of |
---|
150 | $\bar{e}$ through vertical shear, its destruction through stratification, its vertical |
---|
151 | diffusion, and its dissipation of \citet{Kolmogorov1942} type: |
---|
152 | \begin{equation} \label{Eq_zdftke_e} |
---|
153 | \frac{\partial \bar{e}}{\partial t} = |
---|
154 | \frac{K_m}{{e_3}^2 }\;\left[ {\left( {\frac{\partial u}{\partial k}} \right)^2 |
---|
155 | +\left( {\frac{\partial v}{\partial k}} \right)^2} \right] |
---|
156 | -K_\rho\,N^2 |
---|
157 | +\frac{1}{e_3} \;\frac{\partial }{\partial k}\left[ {\frac{A^{vm}}{e_3 } |
---|
158 | \;\frac{\partial \bar{e}}{\partial k}} \right] |
---|
159 | - c_\epsilon \;\frac{\bar {e}^{3/2}}{l_\epsilon } |
---|
160 | \end{equation} |
---|
161 | \begin{equation} \label{Eq_zdftke_kz} |
---|
162 | \begin{split} |
---|
163 | K_m &= C_k\ l_k\ \sqrt {\bar{e}\; } \\ |
---|
164 | K_\rho &= A^{vm} / P_{rt} |
---|
165 | \end{split} |
---|
166 | \end{equation} |
---|
167 | where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}), |
---|
168 | $l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing length scales, |
---|
169 | $P_{rt}$ is the Prandtl number, $K_m$ and $K_\rho$ are the vertical eddy viscosity |
---|
170 | and diffusivity coefficients. The constants $C_k = 0.1$ and $C_\epsilon = \sqrt {2} /2$ |
---|
171 | $\approx 0.7$ are designed to deal with vertical mixing at any depth \citep{Gaspar1990}. |
---|
172 | They are set through namelist parameters \np{nn\_ediff} and \np{nn\_ediss}. |
---|
173 | $P_{rt}$ can be set to unity or, following \citet{Blanke1993}, be a function |
---|
174 | of the local Richardson number, $R_i$: |
---|
175 | \begin{align*} \label{Eq_prt} |
---|
176 | P_{rt} = \begin{cases} |
---|
177 | \ \ \ 1 & \text{if $\ R_i \leq 0.2$} \\ |
---|
178 | 5\,R_i & \text{if $\ 0.2 \leq R_i \leq 2$} \\ |
---|
179 | \ \ 10 & \text{if $\ 2 \leq R_i$} |
---|
180 | \end{cases} |
---|
181 | \end{align*} |
---|
182 | Options are defined through the \ngn{namzdfy\_tke} namelist variables. |
---|
183 | The choice of $P_{rt}$ is controlled by the \np{nn\_pdl} namelist variable. |
---|
184 | |
---|
185 | At the sea surface, the value of $\bar{e}$ is prescribed from the wind |
---|
186 | stress field as $\bar{e}_o = e_{bb} |\tau| / \rho_o$, with $e_{bb}$ the \np{rn\_ebb} |
---|
187 | namelist parameter. The default value of $e_{bb}$ is 3.75. \citep{Gaspar1990}), |
---|
188 | however a much larger value can be used when taking into account the |
---|
189 | surface wave breaking (see below Eq. \eqref{ZDF_Esbc}). |
---|
190 | The bottom value of TKE is assumed to be equal to the value of the level just above. |
---|
191 | The time integration of the $\bar{e}$ equation may formally lead to negative values |
---|
192 | because the numerical scheme does not ensure its positivity. To overcome this |
---|
193 | problem, a cut-off in the minimum value of $\bar{e}$ is used (\np{rn\_emin} |
---|
194 | namelist parameter). Following \citet{Gaspar1990}, the cut-off value is set |
---|
195 | to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$. This allows the subsequent formulations |
---|
196 | to match that of \citet{Gargett1984} for the diffusion in the thermocline and |
---|
197 | deep ocean : $K_\rho = 10^{-3} / N$. |
---|
198 | In addition, a cut-off is applied on $K_m$ and $K_\rho$ to avoid numerical |
---|
199 | instabilities associated with too weak vertical diffusion. They must be |
---|
200 | specified at least larger than the molecular values, and are set through |
---|
201 | \np{rn\_avm0} and \np{rn\_avt0} (namzdf namelist, see \S\ref{ZDF_cst}). |
---|
202 | |
---|
203 | \subsubsection{Turbulent length scale} |
---|
204 | For computational efficiency, the original formulation of the turbulent length |
---|
205 | scales proposed by \citet{Gaspar1990} has been simplified. Four formulations |
---|
206 | are proposed, the choice of which is controlled by the \np{nn\_mxl} namelist |
---|
207 | parameter. The first two are based on the following first order approximation |
---|
208 | \citep{Blanke1993}: |
---|
209 | \begin{equation} \label{Eq_tke_mxl0_1} |
---|
210 | l_k = l_\epsilon = \sqrt {2 \bar{e}\; } / N |
---|
211 | \end{equation} |
---|
212 | which is valid in a stable stratified region with constant values of the Brunt- |
---|
213 | Vais\"{a}l\"{a} frequency. The resulting length scale is bounded by the distance |
---|
214 | to the surface or to the bottom (\np{nn\_mxl} = 0) or by the local vertical scale factor |
---|
215 | (\np{nn\_mxl} = 1). \citet{Blanke1993} notice that this simplification has two major |
---|
216 | drawbacks: it makes no sense for locally unstable stratification and the |
---|
217 | computation no longer uses all the information contained in the vertical density |
---|
218 | profile. To overcome these drawbacks, \citet{Madec1998} introduces the |
---|
219 | \np{nn\_mxl} = 2 or 3 cases, which add an extra assumption concerning the vertical |
---|
220 | gradient of the computed length scale. So, the length scales are first evaluated |
---|
221 | as in \eqref{Eq_tke_mxl0_1} and then bounded such that: |
---|
222 | \begin{equation} \label{Eq_tke_mxl_constraint} |
---|
223 | \frac{1}{e_3 }\left| {\frac{\partial l}{\partial k}} \right| \leq 1 |
---|
224 | \qquad \text{with }\ l = l_k = l_\epsilon |
---|
225 | \end{equation} |
---|
226 | \eqref{Eq_tke_mxl_constraint} means that the vertical variations of the length |
---|
227 | scale cannot be larger than the variations of depth. It provides a better |
---|
228 | approximation of the \citet{Gaspar1990} formulation while being much less |
---|
229 | time consuming. In particular, it allows the length scale to be limited not only |
---|
230 | by the distance to the surface or to the ocean bottom but also by the distance |
---|
231 | to a strongly stratified portion of the water column such as the thermocline |
---|
232 | (Fig.~\ref{Fig_mixing_length}). In order to impose the \eqref{Eq_tke_mxl_constraint} |
---|
233 | constraint, we introduce two additional length scales: $l_{up}$ and $l_{dwn}$, |
---|
234 | the upward and downward length scales, and evaluate the dissipation and |
---|
235 | mixing length scales as (and note that here we use numerical indexing): |
---|
236 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
237 | \begin{figure}[!t] \begin{center} |
---|
238 | \includegraphics[width=1.00\textwidth]{Fig_mixing_length} |
---|
239 | \caption{ \protect\label{Fig_mixing_length} |
---|
240 | Illustration of the mixing length computation. } |
---|
241 | \end{center} |
---|
242 | \end{figure} |
---|
243 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
244 | \begin{equation} \label{Eq_tke_mxl2} |
---|
245 | \begin{aligned} |
---|
246 | l_{up\ \ }^{(k)} &= \min \left( l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3t}^{(k)}\ \ \ \; \right) |
---|
247 | \quad &\text{ from $k=1$ to $jpk$ }\ \\ |
---|
248 | l_{dwn}^{(k)} &= \min \left( l^{(k)} \ , \ l_{dwn}^{(k-1)} + e_{3t}^{(k-1)} \right) |
---|
249 | \quad &\text{ from $k=jpk$ to $1$ }\ \\ |
---|
250 | \end{aligned} |
---|
251 | \end{equation} |
---|
252 | where $l^{(k)}$ is computed using \eqref{Eq_tke_mxl0_1}, |
---|
253 | $i.e.$ $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$. |
---|
254 | |
---|
255 | In the \np{nn\_mxl}~=~2 case, the dissipation and mixing length scales take the same |
---|
256 | value: $ l_k= l_\epsilon = \min \left(\ l_{up} \;,\; l_{dwn}\ \right)$, while in the |
---|
257 | \np{nn\_mxl}~=~3 case, the dissipation and mixing turbulent length scales are give |
---|
258 | as in \citet{Gaspar1990}: |
---|
259 | \begin{equation} \label{Eq_tke_mxl_gaspar} |
---|
260 | \begin{aligned} |
---|
261 | & l_k = \sqrt{\ l_{up} \ \ l_{dwn}\ } \\ |
---|
262 | & l_\epsilon = \min \left(\ l_{up} \;,\; l_{dwn}\ \right) |
---|
263 | \end{aligned} |
---|
264 | \end{equation} |
---|
265 | |
---|
266 | At the ocean surface, a non zero length scale is set through the \np{rn\_mxl0} namelist |
---|
267 | parameter. Usually the surface scale is given by $l_o = \kappa \,z_o$ |
---|
268 | where $\kappa = 0.4$ is von Karman's constant and $z_o$ the roughness |
---|
269 | parameter of the surface. Assuming $z_o=0.1$~m \citep{Craig_Banner_JPO94} |
---|
270 | leads to a 0.04~m, the default value of \np{rn\_mxl0}. In the ocean interior |
---|
271 | a minimum length scale is set to recover the molecular viscosity when $\bar{e}$ |
---|
272 | reach its minimum value ($1.10^{-6}= C_k\, l_{min} \,\sqrt{\bar{e}_{min}}$ ). |
---|
273 | |
---|
274 | |
---|
275 | \subsubsection{Surface wave breaking parameterization} |
---|
276 | %-----------------------------------------------------------------------% |
---|
277 | Following \citet{Mellor_Blumberg_JPO04}, the TKE turbulence closure model has been modified |
---|
278 | to include the effect of surface wave breaking energetics. This results in a reduction of summertime |
---|
279 | surface temperature when the mixed layer is relatively shallow. The \citet{Mellor_Blumberg_JPO04} |
---|
280 | modifications acts on surface length scale and TKE values and air-sea drag coefficient. |
---|
281 | The latter concerns the bulk formulea and is not discussed here. |
---|
282 | |
---|
283 | Following \citet{Craig_Banner_JPO94}, the boundary condition on surface TKE value is : |
---|
284 | \begin{equation} \label{ZDF_Esbc} |
---|
285 | \bar{e}_o = \frac{1}{2}\,\left( 15.8\,\alpha_{CB} \right)^{2/3} \,\frac{|\tau|}{\rho_o} |
---|
286 | \end{equation} |
---|
287 | where $\alpha_{CB}$ is the \citet{Craig_Banner_JPO94} constant of proportionality |
---|
288 | which depends on the ''wave age'', ranging from 57 for mature waves to 146 for |
---|
289 | younger waves \citep{Mellor_Blumberg_JPO04}. |
---|
290 | The boundary condition on the turbulent length scale follows the Charnock's relation: |
---|
291 | \begin{equation} \label{ZDF_Lsbc} |
---|
292 | l_o = \kappa \beta \,\frac{|\tau|}{g\,\rho_o} |
---|
293 | \end{equation} |
---|
294 | where $\kappa=0.40$ is the von Karman constant, and $\beta$ is the Charnock's constant. |
---|
295 | \citet{Mellor_Blumberg_JPO04} suggest $\beta = 2.10^{5}$ the value chosen by \citet{Stacey_JPO99} |
---|
296 | citing observation evidence, and $\alpha_{CB} = 100$ the Craig and Banner's value. |
---|
297 | As the surface boundary condition on TKE is prescribed through $\bar{e}_o = e_{bb} |\tau| / \rho_o$, |
---|
298 | with $e_{bb}$ the \np{rn\_ebb} namelist parameter, setting \np{rn\_ebb}~=~67.83 corresponds |
---|
299 | to $\alpha_{CB} = 100$. Further setting \np{ln\_mxl0} to true applies \eqref{ZDF_Lsbc} |
---|
300 | as surface boundary condition on length scale, with $\beta$ hard coded to the Stacey's value. |
---|
301 | Note that a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) |
---|
302 | is applied on surface $\bar{e}$ value. |
---|
303 | |
---|
304 | |
---|
305 | \subsubsection{Langmuir cells} |
---|
306 | %--------------------------------------% |
---|
307 | Langmuir circulations (LC) can be described as ordered large-scale vertical motions |
---|
308 | in the surface layer of the oceans. Although LC have nothing to do with convection, |
---|
309 | the circulation pattern is rather similar to so-called convective rolls in the atmospheric |
---|
310 | boundary layer. The detailed physics behind LC is described in, for example, |
---|
311 | \citet{Craik_Leibovich_JFM76}. The prevailing explanation is that LC arise from |
---|
312 | a nonlinear interaction between the Stokes drift and wind drift currents. |
---|
313 | |
---|
314 | Here we introduced in the TKE turbulent closure the simple parameterization of |
---|
315 | Langmuir circulations proposed by \citep{Axell_JGR02} for a $k-\epsilon$ turbulent closure. |
---|
316 | The parameterization, tuned against large-eddy simulation, includes the whole effect |
---|
317 | of LC in an extra source terms of TKE, $P_{LC}$. |
---|
318 | The presence of $P_{LC}$ in \eqref{Eq_zdftke_e}, the TKE equation, is controlled |
---|
319 | by setting \np{ln\_lc} to \textit{true} in the namtke namelist. |
---|
320 | |
---|
321 | By making an analogy with the characteristic convective velocity scale |
---|
322 | ($e.g.$, \citet{D'Alessio_al_JPO98}), $P_{LC}$ is assumed to be : |
---|
323 | \begin{equation} |
---|
324 | P_{LC}(z) = \frac{w_{LC}^3(z)}{H_{LC}} |
---|
325 | \end{equation} |
---|
326 | where $w_{LC}(z)$ is the vertical velocity profile of LC, and $H_{LC}$ is the LC depth. |
---|
327 | With no information about the wave field, $w_{LC}$ is assumed to be proportional to |
---|
328 | the Stokes drift $u_s = 0.377\,\,|\tau|^{1/2}$, where $|\tau|$ is the surface wind stress module |
---|
329 | \footnote{Following \citet{Li_Garrett_JMR93}, the surface Stoke drift velocity |
---|
330 | may be expressed as $u_s = 0.016 \,|U_{10m}|$. Assuming an air density of |
---|
331 | $\rho_a=1.22 \,Kg/m^3$ and a drag coefficient of $1.5~10^{-3}$ give the expression |
---|
332 | used of $u_s$ as a function of the module of surface stress}. |
---|
333 | For the vertical variation, $w_{LC}$ is assumed to be zero at the surface as well as |
---|
334 | at a finite depth $H_{LC}$ (which is often close to the mixed layer depth), and simply |
---|
335 | varies as a sine function in between (a first-order profile for the Langmuir cell structures). |
---|
336 | The resulting expression for $w_{LC}$ is : |
---|
337 | \begin{equation} |
---|
338 | w_{LC} = \begin{cases} |
---|
339 | c_{LC} \,u_s \,\sin(- \pi\,z / H_{LC} ) & \text{if $-z \leq H_{LC}$} \\ |
---|
340 | 0 & \text{otherwise} |
---|
341 | \end{cases} |
---|
342 | \end{equation} |
---|
343 | where $c_{LC} = 0.15$ has been chosen by \citep{Axell_JGR02} as a good compromise |
---|
344 | to fit LES data. The chosen value yields maximum vertical velocities $w_{LC}$ of the order |
---|
345 | of a few centimeters per second. The value of $c_{LC}$ is set through the \np{rn\_lc} |
---|
346 | namelist parameter, having in mind that it should stay between 0.15 and 0.54 \citep{Axell_JGR02}. |
---|
347 | |
---|
348 | The $H_{LC}$ is estimated in a similar way as the turbulent length scale of TKE equations: |
---|
349 | $H_{LC}$ is depth to which a water parcel with kinetic energy due to Stoke drift |
---|
350 | can reach on its own by converting its kinetic energy to potential energy, according to |
---|
351 | \begin{equation} |
---|
352 | - \int_{-H_{LC}}^0 { N^2\;z \;dz} = \frac{1}{2} u_s^2 |
---|
353 | \end{equation} |
---|
354 | |
---|
355 | |
---|
356 | \subsubsection{Mixing just below the mixed layer} |
---|
357 | %--------------------------------------------------------------% |
---|
358 | |
---|
359 | Vertical mixing parameterizations commonly used in ocean general circulation models |
---|
360 | tend to produce mixed-layer depths that are too shallow during summer months and windy conditions. |
---|
361 | This bias is particularly acute over the Southern Ocean. |
---|
362 | To overcome this systematic bias, an ad hoc parameterization is introduced into the TKE scheme \cite{Rodgers_2014}. |
---|
363 | The parameterization is an empirical one, $i.e.$ not derived from theoretical considerations, |
---|
364 | but rather is meant to account for observed processes that affect the density structure of |
---|
365 | the ocean’s planetary boundary layer that are not explicitly captured by default in the TKE scheme |
---|
366 | ($i.e.$ near-inertial oscillations and ocean swells and waves). |
---|
367 | |
---|
368 | When using this parameterization ($i.e.$ when \np{nn\_etau}~=~1), the TKE input to the ocean ($S$) |
---|
369 | imposed by the winds in the form of near-inertial oscillations, swell and waves is parameterized |
---|
370 | by \eqref{ZDF_Esbc} the standard TKE surface boundary condition, plus a depth depend one given by: |
---|
371 | \begin{equation} \label{ZDF_Ehtau} |
---|
372 | S = (1-f_i) \; f_r \; e_s \; e^{-z / h_\tau} |
---|
373 | \end{equation} |
---|
374 | where |
---|
375 | $z$ is the depth, |
---|
376 | $e_s$ is TKE surface boundary condition, |
---|
377 | $f_r$ is the fraction of the surface TKE that penetrate in the ocean, |
---|
378 | $h_\tau$ is a vertical mixing length scale that controls exponential shape of the penetration, |
---|
379 | and $f_i$ is the ice concentration (no penetration if $f_i=1$, that is if the ocean is entirely |
---|
380 | covered by sea-ice). |
---|
381 | The value of $f_r$, usually a few percents, is specified through \np{rn\_efr} namelist parameter. |
---|
382 | The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np{nn\_etau}~=~0) |
---|
383 | or a latitude dependent value (varying from 0.5~m at the Equator to a maximum value of 30~m |
---|
384 | at high latitudes (\np{nn\_etau}~=~1). |
---|
385 | |
---|
386 | Note that two other option existe, \np{nn\_etau}~=~2, or 3. They correspond to applying |
---|
387 | \eqref{ZDF_Ehtau} only at the base of the mixed layer, or to using the high frequency part |
---|
388 | of the stress to evaluate the fraction of TKE that penetrate the ocean. |
---|
389 | Those two options are obsolescent features introduced for test purposes. |
---|
390 | They will be removed in the next release. |
---|
391 | |
---|
392 | |
---|
393 | |
---|
394 | % from Burchard et al OM 2008 : |
---|
395 | % the most critical process not reproduced by statistical turbulence models is the activity of |
---|
396 | % internal waves and their interaction with turbulence. After the Reynolds decomposition, |
---|
397 | % internal waves are in principle included in the RANS equations, but later partially |
---|
398 | % excluded by the hydrostatic assumption and the model resolution. |
---|
399 | % Thus far, the representation of internal wave mixing in ocean models has been relatively crude |
---|
400 | % (e.g. Mellor, 1989; Large et al., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002). |
---|
401 | |
---|
402 | |
---|
403 | |
---|
404 | % ------------------------------------------------------------------------------------------------------------- |
---|
405 | % TKE discretization considerations |
---|
406 | % ------------------------------------------------------------------------------------------------------------- |
---|
407 | \subsection{TKE discretization considerations (\protect\key{zdftke})} |
---|
408 | \label{ZDF_tke_ene} |
---|
409 | |
---|
410 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
411 | \begin{figure}[!t] \begin{center} |
---|
412 | \includegraphics[width=1.00\textwidth]{Fig_ZDF_TKE_time_scheme} |
---|
413 | \caption{ \protect\label{Fig_TKE_time_scheme} |
---|
414 | Illustration of the TKE time integration and its links to the momentum and tracer time integration. } |
---|
415 | \end{center} |
---|
416 | \end{figure} |
---|
417 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
418 | |
---|
419 | The production of turbulence by vertical shear (the first term of the right hand side |
---|
420 | of \eqref{Eq_zdftke_e}) should balance the loss of kinetic energy associated with |
---|
421 | the vertical momentum diffusion (first line in \eqref{Eq_PE_zdf}). To do so a special care |
---|
422 | have to be taken for both the time and space discretization of the TKE equation |
---|
423 | \citep{Burchard_OM02,Marsaleix_al_OM08}. |
---|
424 | |
---|
425 | Let us first address the time stepping issue. Fig.~\ref{Fig_TKE_time_scheme} shows |
---|
426 | how the two-level Leap-Frog time stepping of the momentum and tracer equations interplays |
---|
427 | with the one-level forward time stepping of TKE equation. With this framework, the total loss |
---|
428 | of kinetic energy (in 1D for the demonstration) due to the vertical momentum diffusion is |
---|
429 | obtained by multiplying this quantity by $u^t$ and summing the result vertically: |
---|
430 | \begin{equation} \label{Eq_energ1} |
---|
431 | \begin{split} |
---|
432 | \int_{-H}^{\eta} u^t \,\partial_z &\left( {K_m}^t \,(\partial_z u)^{t+\rdt} \right) \,dz \\ |
---|
433 | &= \Bigl[ u^t \,{K_m}^t \,(\partial_z u)^{t+\rdt} \Bigr]_{-H}^{\eta} |
---|
434 | - \int_{-H}^{\eta}{ {K_m}^t \,\partial_z{u^t} \,\partial_z u^{t+\rdt} \,dz } |
---|
435 | \end{split} |
---|
436 | \end{equation} |
---|
437 | Here, the vertical diffusion of momentum is discretized backward in time |
---|
438 | with a coefficient, $K_m$, known at time $t$ (Fig.~\ref{Fig_TKE_time_scheme}), |
---|
439 | as it is required when using the TKE scheme (see \S\ref{STP_forward_imp}). |
---|
440 | The first term of the right hand side of \eqref{Eq_energ1} represents the kinetic energy |
---|
441 | transfer at the surface (atmospheric forcing) and at the bottom (friction effect). |
---|
442 | The second term is always negative. It is the dissipation rate of kinetic energy, |
---|
443 | and thus minus the shear production rate of $\bar{e}$. \eqref{Eq_energ1} |
---|
444 | implies that, to be energetically consistent, the production rate of $\bar{e}$ |
---|
445 | used to compute $(\bar{e})^t$ (and thus ${K_m}^t$) should be expressed as |
---|
446 | ${K_m}^{t-\rdt}\,(\partial_z u)^{t-\rdt} \,(\partial_z u)^t$ (and not by the more straightforward |
---|
447 | $K_m \left( \partial_z u \right)^2$ expression taken at time $t$ or $t-\rdt$). |
---|
448 | |
---|
449 | A similar consideration applies on the destruction rate of $\bar{e}$ due to stratification |
---|
450 | (second term of the right hand side of \eqref{Eq_zdftke_e}). This term |
---|
451 | must balance the input of potential energy resulting from vertical mixing. |
---|
452 | The rate of change of potential energy (in 1D for the demonstration) due vertical |
---|
453 | mixing is obtained by multiplying vertical density diffusion |
---|
454 | tendency by $g\,z$ and and summing the result vertically: |
---|
455 | \begin{equation} \label{Eq_energ2} |
---|
456 | \begin{split} |
---|
457 | \int_{-H}^{\eta} g\,z\,\partial_z &\left( {K_\rho}^t \,(\partial_k \rho)^{t+\rdt} \right) \,dz \\ |
---|
458 | &= \Bigl[ g\,z \,{K_\rho}^t \,(\partial_z \rho)^{t+\rdt} \Bigr]_{-H}^{\eta} |
---|
459 | - \int_{-H}^{\eta}{ g \,{K_\rho}^t \,(\partial_k \rho)^{t+\rdt} } \,dz \\ |
---|
460 | &= - \Bigl[ z\,{K_\rho}^t \,(N^2)^{t+\rdt} \Bigr]_{-H}^{\eta} |
---|
461 | + \int_{-H}^{\eta}{ \rho^{t+\rdt} \, {K_\rho}^t \,(N^2)^{t+\rdt} \,dz } |
---|
462 | \end{split} |
---|
463 | \end{equation} |
---|
464 | where we use $N^2 = -g \,\partial_k \rho / (e_3 \rho)$. |
---|
465 | The first term of the right hand side of \eqref{Eq_energ2} is always zero |
---|
466 | because there is no diffusive flux through the ocean surface and bottom). |
---|
467 | The second term is minus the destruction rate of $\bar{e}$ due to stratification. |
---|
468 | Therefore \eqref{Eq_energ1} implies that, to be energetically consistent, the product |
---|
469 | ${K_\rho}^{t-\rdt}\,(N^2)^t$ should be used in \eqref{Eq_zdftke_e}, the TKE equation. |
---|
470 | |
---|
471 | Let us now address the space discretization issue. |
---|
472 | The vertical eddy coefficients are defined at $w$-point whereas the horizontal velocity |
---|
473 | components are in the centre of the side faces of a $t$-box in staggered C-grid |
---|
474 | (Fig.\ref{Fig_cell}). A space averaging is thus required to obtain the shear TKE production term. |
---|
475 | By redoing the \eqref{Eq_energ1} in the 3D case, it can be shown that the product of |
---|
476 | eddy coefficient by the shear at $t$ and $t-\rdt$ must be performed prior to the averaging. |
---|
477 | Furthermore, the possible time variation of $e_3$ (\key{vvl} case) have to be taken into |
---|
478 | account. |
---|
479 | |
---|
480 | The above energetic considerations leads to |
---|
481 | the following final discrete form for the TKE equation: |
---|
482 | \begin{equation} \label{Eq_zdftke_ene} |
---|
483 | \begin{split} |
---|
484 | \frac { (\bar{e})^t - (\bar{e})^{t-\rdt} } {\rdt} \equiv |
---|
485 | \Biggl\{ \Biggr. |
---|
486 | &\overline{ \left( \left(\overline{K_m}^{\,i+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[u^{t+\rdt}]}{{e_3u}^{t+\rdt} } |
---|
487 | \ \frac{\delta_{k+1/2}[u^ t ]}{{e_3u}^ t } \right) }^{\,i} \\ |
---|
488 | +&\overline{ \left( \left(\overline{K_m}^{\,j+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[v^{t+\rdt}]}{{e_3v}^{t+\rdt} } |
---|
489 | \ \frac{\delta_{k+1/2}[v^ t ]}{{e_3v}^ t } \right) }^{\,j} |
---|
490 | \Biggr. \Biggr\} \\ |
---|
491 | % |
---|
492 | - &{K_\rho}^{t-\rdt}\,{(N^2)^t} \\ |
---|
493 | % |
---|
494 | +&\frac{1}{{e_3w}^{t+\rdt}} \;\delta_{k+1/2} \left[ {K_m}^{t-\rdt} \,\frac{\delta_{k}[(\bar{e})^{t+\rdt}]} {{e_3w}^{t+\rdt}} \right] \\ |
---|
495 | % |
---|
496 | - &c_\epsilon \; \left( \frac{\sqrt{\bar {e}}}{l_\epsilon}\right)^{t-\rdt}\,(\bar {e})^{t+\rdt} |
---|
497 | \end{split} |
---|
498 | \end{equation} |
---|
499 | where the last two terms in \eqref{Eq_zdftke_ene} (vertical diffusion and Kolmogorov dissipation) |
---|
500 | are time stepped using a backward scheme (see\S\ref{STP_forward_imp}). |
---|
501 | Note that the Kolmogorov term has been linearized in time in order to render |
---|
502 | the implicit computation possible. The restart of the TKE scheme |
---|
503 | requires the storage of $\bar {e}$, $K_m$, $K_\rho$ and $l_\epsilon$ as they all appear in |
---|
504 | the right hand side of \eqref{Eq_zdftke_ene}. For the latter, it is in fact |
---|
505 | the ratio $\sqrt{\bar{e}}/l_\epsilon$ which is stored. |
---|
506 | |
---|
507 | % ------------------------------------------------------------------------------------------------------------- |
---|
508 | % GLS Generic Length Scale Scheme |
---|
509 | % ------------------------------------------------------------------------------------------------------------- |
---|
510 | \subsection{GLS Generic Length Scale (\protect\key{zdfgls})} |
---|
511 | \label{ZDF_gls} |
---|
512 | |
---|
513 | %--------------------------------------------namzdf_gls--------------------------------------------------------- |
---|
514 | \namdisplay{namzdf_gls} |
---|
515 | %-------------------------------------------------------------------------------------------------------------- |
---|
516 | |
---|
517 | The Generic Length Scale (GLS) scheme is a turbulent closure scheme based on |
---|
518 | two prognostic equations: one for the turbulent kinetic energy $\bar {e}$, and another |
---|
519 | for the generic length scale, $\psi$ \citep{Umlauf_Burchard_JMS03, Umlauf_Burchard_CSR05}. |
---|
520 | This later variable is defined as : $\psi = {C_{0\mu}}^{p} \ {\bar{e}}^{m} \ l^{n}$, |
---|
521 | where the triplet $(p, m, n)$ value given in Tab.\ref{Tab_GLS} allows to recover |
---|
522 | a number of well-known turbulent closures ($k$-$kl$ \citep{Mellor_Yamada_1982}, |
---|
523 | $k$-$\epsilon$ \citep{Rodi_1987}, $k$-$\omega$ \citep{Wilcox_1988} |
---|
524 | among others \citep{Umlauf_Burchard_JMS03,Kantha_Carniel_CSR05}). |
---|
525 | The GLS scheme is given by the following set of equations: |
---|
526 | \begin{equation} \label{Eq_zdfgls_e} |
---|
527 | \frac{\partial \bar{e}}{\partial t} = |
---|
528 | \frac{K_m}{\sigma_e e_3 }\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2 |
---|
529 | +\left( \frac{\partial v}{\partial k} \right)^2} \right] |
---|
530 | -K_\rho \,N^2 |
---|
531 | +\frac{1}{e_3}\,\frac{\partial}{\partial k} \left[ \frac{K_m}{e_3}\,\frac{\partial \bar{e}}{\partial k} \right] |
---|
532 | - \epsilon |
---|
533 | \end{equation} |
---|
534 | |
---|
535 | \begin{equation} \label{Eq_zdfgls_psi} |
---|
536 | \begin{split} |
---|
537 | \frac{\partial \psi}{\partial t} =& \frac{\psi}{\bar{e}} \left\{ |
---|
538 | \frac{C_1\,K_m}{\sigma_{\psi} {e_3}}\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2 |
---|
539 | +\left( \frac{\partial v}{\partial k} \right)^2} \right] |
---|
540 | - C_3 \,K_\rho\,N^2 - C_2 \,\epsilon \,Fw \right\} \\ |
---|
541 | &+\frac{1}{e_3} \;\frac{\partial }{\partial k}\left[ {\frac{K_m}{e_3 } |
---|
542 | \;\frac{\partial \psi}{\partial k}} \right]\; |
---|
543 | \end{split} |
---|
544 | \end{equation} |
---|
545 | |
---|
546 | \begin{equation} \label{Eq_zdfgls_kz} |
---|
547 | \begin{split} |
---|
548 | K_m &= C_{\mu} \ \sqrt {\bar{e}} \ l \\ |
---|
549 | K_\rho &= C_{\mu'}\ \sqrt {\bar{e}} \ l |
---|
550 | \end{split} |
---|
551 | \end{equation} |
---|
552 | |
---|
553 | \begin{equation} \label{Eq_zdfgls_eps} |
---|
554 | {\epsilon} = C_{0\mu} \,\frac{\bar {e}^{3/2}}{l} \; |
---|
555 | \end{equation} |
---|
556 | where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}) |
---|
557 | and $\epsilon$ the dissipation rate. |
---|
558 | The constants $C_1$, $C_2$, $C_3$, ${\sigma_e}$, ${\sigma_{\psi}}$ and the wall function ($Fw$) |
---|
559 | depends of the choice of the turbulence model. Four different turbulent models are pre-defined |
---|
560 | (Tab.\ref{Tab_GLS}). They are made available through the \np{nn\_clo} namelist parameter. |
---|
561 | |
---|
562 | %--------------------------------------------------TABLE-------------------------------------------------- |
---|
563 | \begin{table}[htbp] \begin{center} |
---|
564 | %\begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}c} |
---|
565 | \begin{tabular}{ccccc} |
---|
566 | & $k-kl$ & $k-\epsilon$ & $k-\omega$ & generic \\ |
---|
567 | % & \citep{Mellor_Yamada_1982} & \citep{Rodi_1987} & \citep{Wilcox_1988} & \\ |
---|
568 | \hline \hline |
---|
569 | \np{nn\_clo} & \textbf{0} & \textbf{1} & \textbf{2} & \textbf{3} \\ |
---|
570 | \hline |
---|
571 | $( p , n , m )$ & ( 0 , 1 , 1 ) & ( 3 , 1.5 , -1 ) & ( -1 , 0.5 , -1 ) & ( 2 , 1 , -0.67 ) \\ |
---|
572 | $\sigma_k$ & 2.44 & 1. & 2. & 0.8 \\ |
---|
573 | $\sigma_\psi$ & 2.44 & 1.3 & 2. & 1.07 \\ |
---|
574 | $C_1$ & 0.9 & 1.44 & 0.555 & 1. \\ |
---|
575 | $C_2$ & 0.5 & 1.92 & 0.833 & 1.22 \\ |
---|
576 | $C_3$ & 1. & 1. & 1. & 1. \\ |
---|
577 | $F_{wall}$ & Yes & -- & -- & -- \\ |
---|
578 | \hline |
---|
579 | \hline |
---|
580 | \end{tabular} |
---|
581 | \caption{ \protect\label{Tab_GLS} |
---|
582 | Set of predefined GLS parameters, or equivalently predefined turbulence models available |
---|
583 | with \protect\key{zdfgls} and controlled by the \protect\np{nn\_clos} namelist variable in \protect\ngn{namzdf\_gls} .} |
---|
584 | \end{center} \end{table} |
---|
585 | %-------------------------------------------------------------------------------------------------------------- |
---|
586 | |
---|
587 | In the Mellor-Yamada model, the negativity of $n$ allows to use a wall function to force |
---|
588 | the convergence of the mixing length towards $K z_b$ ($K$: Kappa and $z_b$: rugosity length) |
---|
589 | value near physical boundaries (logarithmic boundary layer law). $C_{\mu}$ and $C_{\mu'}$ |
---|
590 | are calculated from stability function proposed by \citet{Galperin_al_JAS88}, or by \citet{Kantha_Clayson_1994} |
---|
591 | or one of the two functions suggested by \citet{Canuto_2001} (\np{nn\_stab\_func} = 0, 1, 2 or 3, resp.). |
---|
592 | The value of $C_{0\mu}$ depends of the choice of the stability function. |
---|
593 | |
---|
594 | The surface and bottom boundary condition on both $\bar{e}$ and $\psi$ can be calculated |
---|
595 | thanks to Dirichlet or Neumann condition through \np{nn\_tkebc\_surf} and \np{nn\_tkebc\_bot}, resp. |
---|
596 | As for TKE closure , the wave effect on the mixing is considered when \np{ln\_crban}~=~true |
---|
597 | \citep{Craig_Banner_JPO94, Mellor_Blumberg_JPO04}. The \np{rn\_crban} namelist parameter |
---|
598 | is $\alpha_{CB}$ in \eqref{ZDF_Esbc} and \np{rn\_charn} provides the value of $\beta$ in \eqref{ZDF_Lsbc}. |
---|
599 | |
---|
600 | The $\psi$ equation is known to fail in stably stratified flows, and for this reason |
---|
601 | almost all authors apply a clipping of the length scale as an \textit{ad hoc} remedy. |
---|
602 | With this clipping, the maximum permissible length scale is determined by |
---|
603 | $l_{max} = c_{lim} \sqrt{2\bar{e}}/ N$. A value of $c_{lim} = 0.53$ is often used |
---|
604 | \citep{Galperin_al_JAS88}. \cite{Umlauf_Burchard_CSR05} show that the value of |
---|
605 | the clipping factor is of crucial importance for the entrainment depth predicted in |
---|
606 | stably stratified situations, and that its value has to be chosen in accordance |
---|
607 | with the algebraic model for the turbulent fluxes. The clipping is only activated |
---|
608 | if \np{ln\_length\_lim}=true, and the $c_{lim}$ is set to the \np{rn\_clim\_galp} value. |
---|
609 | |
---|
610 | The time and space discretization of the GLS equations follows the same energetic |
---|
611 | consideration as for the TKE case described in \S\ref{ZDF_tke_ene} \citep{Burchard_OM02}. |
---|
612 | Examples of performance of the 4 turbulent closure scheme can be found in \citet{Warner_al_OM05}. |
---|
613 | |
---|
614 | % ------------------------------------------------------------------------------------------------------------- |
---|
615 | % OSM OSMOSIS BL Scheme |
---|
616 | % ------------------------------------------------------------------------------------------------------------- |
---|
617 | \subsection{OSM OSMOSIS Boundary Layer scheme (\protect\key{zdfosm})} |
---|
618 | \label{ZDF_osm} |
---|
619 | |
---|
620 | %--------------------------------------------namzdf_osm--------------------------------------------------------- |
---|
621 | \namdisplay{namzdf_osm} |
---|
622 | %-------------------------------------------------------------------------------------------------------------- |
---|
623 | |
---|
624 | The OSMOSIS turbulent closure scheme is based on...... TBC |
---|
625 | |
---|
626 | % ================================================================ |
---|
627 | % Convection |
---|
628 | % ================================================================ |
---|
629 | \section{Convection} |
---|
630 | \label{ZDF_conv} |
---|
631 | |
---|
632 | %--------------------------------------------namzdf-------------------------------------------------------- |
---|
633 | \namdisplay{namzdf} |
---|
634 | %-------------------------------------------------------------------------------------------------------------- |
---|
635 | |
---|
636 | Static instabilities (i.e. light potential densities under heavy ones) may |
---|
637 | occur at particular ocean grid points. In nature, convective processes |
---|
638 | quickly re-establish the static stability of the water column. These |
---|
639 | processes have been removed from the model via the hydrostatic |
---|
640 | assumption so they must be parameterized. Three parameterisations |
---|
641 | are available to deal with convective processes: a non-penetrative |
---|
642 | convective adjustment or an enhanced vertical diffusion, or/and the |
---|
643 | use of a turbulent closure scheme. |
---|
644 | |
---|
645 | % ------------------------------------------------------------------------------------------------------------- |
---|
646 | % Non-Penetrative Convective Adjustment |
---|
647 | % ------------------------------------------------------------------------------------------------------------- |
---|
648 | \subsection [Non-Penetrative Convective Adjustment (\protect\np{ln\_tranpc}) ] |
---|
649 | {Non-Penetrative Convective Adjustment (\protect\np{ln\_tranpc}=.true.) } |
---|
650 | \label{ZDF_npc} |
---|
651 | |
---|
652 | %--------------------------------------------namzdf-------------------------------------------------------- |
---|
653 | \namdisplay{namzdf} |
---|
654 | %-------------------------------------------------------------------------------------------------------------- |
---|
655 | |
---|
656 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
657 | \begin{figure}[!htb] \begin{center} |
---|
658 | \includegraphics[width=0.90\textwidth]{Fig_npc} |
---|
659 | \caption{ \protect\label{Fig_npc} |
---|
660 | Example of an unstable density profile treated by the non penetrative |
---|
661 | convective adjustment algorithm. $1^{st}$ step: the initial profile is checked from |
---|
662 | the surface to the bottom. It is found to be unstable between levels 3 and 4. |
---|
663 | They are mixed. The resulting $\rho$ is still larger than $\rho$(5): levels 3 to 5 |
---|
664 | are mixed. The resulting $\rho$ is still larger than $\rho$(6): levels 3 to 6 are |
---|
665 | mixed. The $1^{st}$ step ends since the density profile is then stable below |
---|
666 | the level 3. $2^{nd}$ step: the new $\rho$ profile is checked following the same |
---|
667 | procedure as in $1^{st}$ step: levels 2 to 5 are mixed. The new density profile |
---|
668 | is checked. It is found stable: end of algorithm.} |
---|
669 | \end{center} \end{figure} |
---|
670 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
671 | |
---|
672 | Options are defined through the \ngn{namzdf} namelist variables. |
---|
673 | The non-penetrative convective adjustment is used when \np{ln\_zdfnpc}~=~\textit{true}. |
---|
674 | It is applied at each \np{nn\_npc} time step and mixes downwards instantaneously |
---|
675 | the statically unstable portion of the water column, but only until the density |
---|
676 | structure becomes neutrally stable ($i.e.$ until the mixed portion of the water |
---|
677 | column has \textit{exactly} the density of the water just below) \citep{Madec_al_JPO91}. |
---|
678 | The associated algorithm is an iterative process used in the following way |
---|
679 | (Fig. \ref{Fig_npc}): starting from the top of the ocean, the first instability is |
---|
680 | found. Assume in the following that the instability is located between levels |
---|
681 | $k$ and $k+1$. The temperature and salinity in the two levels are |
---|
682 | vertically mixed, conserving the heat and salt contents of the water column. |
---|
683 | The new density is then computed by a linear approximation. If the new |
---|
684 | density profile is still unstable between levels $k+1$ and $k+2$, levels $k$, |
---|
685 | $k+1$ and $k+2$ are then mixed. This process is repeated until stability is |
---|
686 | established below the level $k$ (the mixing process can go down to the |
---|
687 | ocean bottom). The algorithm is repeated to check if the density profile |
---|
688 | between level $k-1$ and $k$ is unstable and/or if there is no deeper instability. |
---|
689 | |
---|
690 | This algorithm is significantly different from mixing statically unstable levels |
---|
691 | two by two. The latter procedure cannot converge with a finite number |
---|
692 | of iterations for some vertical profiles while the algorithm used in \NEMO |
---|
693 | converges for any profile in a number of iterations which is less than the |
---|
694 | number of vertical levels. This property is of paramount importance as |
---|
695 | pointed out by \citet{Killworth1989}: it avoids the existence of permanent |
---|
696 | and unrealistic static instabilities at the sea surface. This non-penetrative |
---|
697 | convective algorithm has been proved successful in studies of the deep |
---|
698 | water formation in the north-western Mediterranean Sea |
---|
699 | \citep{Madec_al_JPO91, Madec_al_DAO91, Madec_Crepon_Bk91}. |
---|
700 | |
---|
701 | The current implementation has been modified in order to deal with any non linear |
---|
702 | equation of seawater (L. Brodeau, personnal communication). |
---|
703 | Two main differences have been introduced compared to the original algorithm: |
---|
704 | $(i)$ the stability is now checked using the Brunt-V\"{a}is\"{a}l\"{a} frequency |
---|
705 | (not the the difference in potential density) ; |
---|
706 | $(ii)$ when two levels are found unstable, their thermal and haline expansion coefficients |
---|
707 | are vertically mixed in the same way their temperature and salinity has been mixed. |
---|
708 | These two modifications allow the algorithm to perform properly and accurately |
---|
709 | with TEOS10 or EOS-80 without having to recompute the expansion coefficients at each |
---|
710 | mixing iteration. |
---|
711 | |
---|
712 | % ------------------------------------------------------------------------------------------------------------- |
---|
713 | % Enhanced Vertical Diffusion |
---|
714 | % ------------------------------------------------------------------------------------------------------------- |
---|
715 | \subsection [Enhanced Vertical Diffusion (\protect\np{ln\_zdfevd})] |
---|
716 | {Enhanced Vertical Diffusion (\protect\np{ln\_zdfevd}=true)} |
---|
717 | \label{ZDF_evd} |
---|
718 | |
---|
719 | %--------------------------------------------namzdf-------------------------------------------------------- |
---|
720 | \namdisplay{namzdf} |
---|
721 | %-------------------------------------------------------------------------------------------------------------- |
---|
722 | |
---|
723 | Options are defined through the \ngn{namzdf} namelist variables. |
---|
724 | The enhanced vertical diffusion parameterisation is used when \np{ln\_zdfevd}=true. |
---|
725 | In this case, the vertical eddy mixing coefficients are assigned very large values |
---|
726 | (a typical value is $10\;m^2s^{-1})$ in regions where the stratification is unstable |
---|
727 | ($i.e.$ when $N^2$ the Brunt-Vais\"{a}l\"{a} frequency is negative) |
---|
728 | \citep{Lazar_PhD97, Lazar_al_JPO99}. This is done either on tracers only |
---|
729 | (\np{nn\_evdm}=0) or on both momentum and tracers (\np{nn\_evdm}=1). |
---|
730 | |
---|
731 | In practice, where $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$, and |
---|
732 | if \np{nn\_evdm}=1, the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm}$ |
---|
733 | values also, are set equal to the namelist parameter \np{rn\_avevd}. A typical value |
---|
734 | for $rn\_avevd$ is between 1 and $100~m^2.s^{-1}$. This parameterisation of |
---|
735 | convective processes is less time consuming than the convective adjustment |
---|
736 | algorithm presented above when mixing both tracers and momentum in the |
---|
737 | case of static instabilities. It requires the use of an implicit time stepping on |
---|
738 | vertical diffusion terms (i.e. \np{ln\_zdfexp}=false). |
---|
739 | |
---|
740 | Note that the stability test is performed on both \textit{before} and \textit{now} |
---|
741 | values of $N^2$. This removes a potential source of divergence of odd and |
---|
742 | even time step in a leapfrog environment \citep{Leclair_PhD2010} (see \S\ref{STP_mLF}). |
---|
743 | |
---|
744 | % ------------------------------------------------------------------------------------------------------------- |
---|
745 | % Turbulent Closure Scheme |
---|
746 | % ------------------------------------------------------------------------------------------------------------- |
---|
747 | \subsection[Turbulent Closure Scheme (\protect\key{zdf\{tke, gls, osm\}})]{Turbulent Closure Scheme (\protect\key{zdftke}, \protect\key{zdfgls} or \protect\key{zdfosm})} |
---|
748 | \label{ZDF_tcs} |
---|
749 | |
---|
750 | The turbulent closure scheme presented in \S\ref{ZDF_tke} and \S\ref{ZDF_gls} |
---|
751 | (\key{zdftke} or \key{zdftke} is defined) in theory solves the problem of statically |
---|
752 | unstable density profiles. In such a case, the term corresponding to the |
---|
753 | destruction of turbulent kinetic energy through stratification in \eqref{Eq_zdftke_e} |
---|
754 | or \eqref{Eq_zdfgls_e} becomes a source term, since $N^2$ is negative. |
---|
755 | It results in large values of $A_T^{vT}$ and $A_T^{vT}$, and also the four neighbouring |
---|
756 | $A_u^{vm} {and}\;A_v^{vm}$ (up to $1\;m^2s^{-1})$. These large values |
---|
757 | restore the static stability of the water column in a way similar to that of the |
---|
758 | enhanced vertical diffusion parameterisation (\S\ref{ZDF_evd}). However, |
---|
759 | in the vicinity of the sea surface (first ocean layer), the eddy coefficients |
---|
760 | computed by the turbulent closure scheme do not usually exceed $10^{-2}m.s^{-1}$, |
---|
761 | because the mixing length scale is bounded by the distance to the sea surface. |
---|
762 | It can thus be useful to combine the enhanced vertical |
---|
763 | diffusion with the turbulent closure scheme, $i.e.$ setting the \np{ln\_zdfnpc} |
---|
764 | namelist parameter to true and defining the turbulent closure CPP key all together. |
---|
765 | |
---|
766 | The KPP turbulent closure scheme already includes enhanced vertical diffusion |
---|
767 | in the case of convection, as governed by the variables $bvsqcon$ and $difcon$ |
---|
768 | found in \mdl{zdfkpp}, therefore \np{ln\_zdfevd}=false should be used with the KPP |
---|
769 | scheme. %gm% + one word on non local flux with KPP scheme trakpp.F90 module... |
---|
770 | |
---|
771 | % ================================================================ |
---|
772 | % Double Diffusion Mixing |
---|
773 | % ================================================================ |
---|
774 | \section [Double Diffusion Mixing (\protect\key{zdfddm})] |
---|
775 | {Double Diffusion Mixing (\protect\key{zdfddm})} |
---|
776 | \label{ZDF_ddm} |
---|
777 | |
---|
778 | %-------------------------------------------namzdf_ddm------------------------------------------------- |
---|
779 | \namdisplay{namzdf_ddm} |
---|
780 | %-------------------------------------------------------------------------------------------------------------- |
---|
781 | |
---|
782 | Options are defined through the \ngn{namzdf\_ddm} namelist variables. |
---|
783 | Double diffusion occurs when relatively warm, salty water overlies cooler, fresher |
---|
784 | water, or vice versa. The former condition leads to salt fingering and the latter |
---|
785 | to diffusive convection. Double-diffusive phenomena contribute to diapycnal |
---|
786 | mixing in extensive regions of the ocean. \citet{Merryfield1999} include a |
---|
787 | parameterisation of such phenomena in a global ocean model and show that |
---|
788 | it leads to relatively minor changes in circulation but exerts significant regional |
---|
789 | influences on temperature and salinity. This parameterisation has been |
---|
790 | introduced in \mdl{zdfddm} module and is controlled by the \key{zdfddm} CPP key. |
---|
791 | |
---|
792 | Diapycnal mixing of S and T are described by diapycnal diffusion coefficients |
---|
793 | \begin{align*} % \label{Eq_zdfddm_Kz} |
---|
794 | &A^{vT} = A_o^{vT}+A_f^{vT}+A_d^{vT} \\ |
---|
795 | &A^{vS} = A_o^{vS}+A_f^{vS}+A_d^{vS} |
---|
796 | \end{align*} |
---|
797 | where subscript $f$ represents mixing by salt fingering, $d$ by diffusive convection, |
---|
798 | and $o$ by processes other than double diffusion. The rates of double-diffusive |
---|
799 | mixing depend on the buoyancy ratio $R_\rho = \alpha \partial_z T / \beta \partial_z S$, |
---|
800 | where $\alpha$ and $\beta$ are coefficients of thermal expansion and saline |
---|
801 | contraction (see \S\ref{TRA_eos}). To represent mixing of $S$ and $T$ by salt |
---|
802 | fingering, we adopt the diapycnal diffusivities suggested by Schmitt (1981): |
---|
803 | \begin{align} \label{Eq_zdfddm_f} |
---|
804 | A_f^{vS} &= \begin{cases} |
---|
805 | \frac{A^{\ast v}}{1+(R_\rho / R_c)^n } &\text{if $R_\rho > 1$ and $N^2>0$ } \\ |
---|
806 | 0 &\text{otherwise} |
---|
807 | \end{cases} |
---|
808 | \\ \label{Eq_zdfddm_f_T} |
---|
809 | A_f^{vT} &= 0.7 \ A_f^{vS} / R_\rho |
---|
810 | \end{align} |
---|
811 | |
---|
812 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
813 | \begin{figure}[!t] \begin{center} |
---|
814 | \includegraphics[width=0.99\textwidth]{Fig_zdfddm} |
---|
815 | \caption{ \protect\label{Fig_zdfddm} |
---|
816 | From \citet{Merryfield1999} : (a) Diapycnal diffusivities $A_f^{vT}$ |
---|
817 | and $A_f^{vS}$ for temperature and salt in regions of salt fingering. Heavy |
---|
818 | curves denote $A^{\ast v} = 10^{-3}~m^2.s^{-1}$ and thin curves |
---|
819 | $A^{\ast v} = 10^{-4}~m^2.s^{-1}$ ; (b) diapycnal diffusivities $A_d^{vT}$ and |
---|
820 | $A_d^{vS}$ for temperature and salt in regions of diffusive convection. Heavy |
---|
821 | curves denote the Federov parameterisation and thin curves the Kelley |
---|
822 | parameterisation. The latter is not implemented in \NEMO. } |
---|
823 | \end{center} \end{figure} |
---|
824 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
825 | |
---|
826 | The factor 0.7 in \eqref{Eq_zdfddm_f_T} reflects the measured ratio |
---|
827 | $\alpha F_T /\beta F_S \approx 0.7$ of buoyancy flux of heat to buoyancy |
---|
828 | flux of salt ($e.g.$, \citet{McDougall_Taylor_JMR84}). Following \citet{Merryfield1999}, |
---|
829 | we adopt $R_c = 1.6$, $n = 6$, and $A^{\ast v} = 10^{-4}~m^2.s^{-1}$. |
---|
830 | |
---|
831 | To represent mixing of S and T by diffusive layering, the diapycnal diffusivities suggested by Federov (1988) is used: |
---|
832 | \begin{align} \label{Eq_zdfddm_d} |
---|
833 | A_d^{vT} &= \begin{cases} |
---|
834 | 1.3635 \, \exp{\left( 4.6\, \exp{ \left[ -0.54\,( R_{\rho}^{-1} - 1 ) \right] } \right)} |
---|
835 | &\text{if $0<R_\rho < 1$ and $N^2>0$ } \\ |
---|
836 | 0 &\text{otherwise} |
---|
837 | \end{cases} |
---|
838 | \\ \label{Eq_zdfddm_d_S} |
---|
839 | A_d^{vS} &= \begin{cases} |
---|
840 | A_d^{vT}\ \left( 1.85\,R_{\rho} - 0.85 \right) |
---|
841 | &\text{if $0.5 \leq R_\rho<1$ and $N^2>0$ } \\ |
---|
842 | A_d^{vT} \ 0.15 \ R_\rho |
---|
843 | &\text{if $\ \ 0 < R_\rho<0.5$ and $N^2>0$ } \\ |
---|
844 | 0 &\text{otherwise} |
---|
845 | \end{cases} |
---|
846 | \end{align} |
---|
847 | |
---|
848 | The dependencies of \eqref{Eq_zdfddm_f} to \eqref{Eq_zdfddm_d_S} on $R_\rho$ |
---|
849 | are illustrated in Fig.~\ref{Fig_zdfddm}. Implementing this requires computing |
---|
850 | $R_\rho$ at each grid point on every time step. This is done in \mdl{eosbn2} at the |
---|
851 | same time as $N^2$ is computed. This avoids duplication in the computation of |
---|
852 | $\alpha$ and $\beta$ (which is usually quite expensive). |
---|
853 | |
---|
854 | % ================================================================ |
---|
855 | % Bottom Friction |
---|
856 | % ================================================================ |
---|
857 | \section [Bottom and Top Friction (\textit{zdfbfr})] {Bottom and Top Friction (\protect\mdl{zdfbfr} module)} |
---|
858 | \label{ZDF_bfr} |
---|
859 | |
---|
860 | %--------------------------------------------nambfr-------------------------------------------------------- |
---|
861 | \namdisplay{nambfr} |
---|
862 | %-------------------------------------------------------------------------------------------------------------- |
---|
863 | |
---|
864 | Options to define the top and bottom friction are defined through the \ngn{nambfr} namelist variables. |
---|
865 | The bottom friction represents the friction generated by the bathymetry. |
---|
866 | The top friction represents the friction generated by the ice shelf/ocean interface. |
---|
867 | As the friction processes at the top and bottom are treated in similar way, |
---|
868 | only the bottom friction is described in detail below. |
---|
869 | |
---|
870 | |
---|
871 | Both the surface momentum flux (wind stress) and the bottom momentum |
---|
872 | flux (bottom friction) enter the equations as a condition on the vertical |
---|
873 | diffusive flux. For the bottom boundary layer, one has: |
---|
874 | \begin{equation} \label{Eq_zdfbfr_flux} |
---|
875 | A^{vm} \left( \partial {\textbf U}_h / \partial z \right) = {{\cal F}}_h^{\textbf U} |
---|
876 | \end{equation} |
---|
877 | where ${\cal F}_h^{\textbf U}$ is represents the downward flux of horizontal momentum |
---|
878 | outside the logarithmic turbulent boundary layer (thickness of the order of |
---|
879 | 1~m in the ocean). How ${\cal F}_h^{\textbf U}$ influences the interior depends on the |
---|
880 | vertical resolution of the model near the bottom relative to the Ekman layer |
---|
881 | depth. For example, in order to obtain an Ekman layer depth |
---|
882 | $d = \sqrt{2\;A^{vm}} / f = 50$~m, one needs a vertical diffusion coefficient |
---|
883 | $A^{vm} = 0.125$~m$^2$s$^{-1}$ (for a Coriolis frequency |
---|
884 | $f = 10^{-4}$~m$^2$s$^{-1}$). With a background diffusion coefficient |
---|
885 | $A^{vm} = 10^{-4}$~m$^2$s$^{-1}$, the Ekman layer depth is only 1.4~m. |
---|
886 | When the vertical mixing coefficient is this small, using a flux condition is |
---|
887 | equivalent to entering the viscous forces (either wind stress or bottom friction) |
---|
888 | as a body force over the depth of the top or bottom model layer. To illustrate |
---|
889 | this, consider the equation for $u$ at $k$, the last ocean level: |
---|
890 | \begin{equation} \label{Eq_zdfbfr_flux2} |
---|
891 | \frac{\partial u_k}{\partial t} = \frac{1}{e_{3u}} \left[ \frac{A_{uw}^{vm}}{e_{3uw}} \delta_{k+1/2}\;[u] - {\cal F}^u_h \right] \approx - \frac{{\cal F}^u_{h}}{e_{3u}} |
---|
892 | \end{equation} |
---|
893 | If the bottom layer thickness is 200~m, the Ekman transport will |
---|
894 | be distributed over that depth. On the other hand, if the vertical resolution |
---|
895 | is high (1~m or less) and a turbulent closure model is used, the turbulent |
---|
896 | Ekman layer will be represented explicitly by the model. However, the |
---|
897 | logarithmic layer is never represented in current primitive equation model |
---|
898 | applications: it is \emph{necessary} to parameterize the flux ${\cal F}^u_h $. |
---|
899 | Two choices are available in \NEMO: a linear and a quadratic bottom friction. |
---|
900 | Note that in both cases, the rotation between the interior velocity and the |
---|
901 | bottom friction is neglected in the present release of \NEMO. |
---|
902 | |
---|
903 | In the code, the bottom friction is imposed by adding the trend due to the bottom |
---|
904 | friction to the general momentum trend in \mdl{dynbfr}. For the time-split surface |
---|
905 | pressure gradient algorithm, the momentum trend due to the barotropic component |
---|
906 | needs to be handled separately. For this purpose it is convenient to compute and |
---|
907 | store coefficients which can be simply combined with bottom velocities and geometric |
---|
908 | values to provide the momentum trend due to bottom friction. |
---|
909 | These coefficients are computed in \mdl{zdfbfr} and generally take the form |
---|
910 | $c_b^{\textbf U}$ where: |
---|
911 | \begin{equation} \label{Eq_zdfbfr_bdef} |
---|
912 | \frac{\partial {\textbf U_h}}{\partial t} = |
---|
913 | - \frac{{\cal F}^{\textbf U}_{h}}{e_{3u}} = \frac{c_b^{\textbf U}}{e_{3u}} \;{\textbf U}_h^b |
---|
914 | \end{equation} |
---|
915 | where $\textbf{U}_h^b = (u_b\;,\;v_b)$ is the near-bottom, horizontal, ocean velocity. |
---|
916 | |
---|
917 | % ------------------------------------------------------------------------------------------------------------- |
---|
918 | % Linear Bottom Friction |
---|
919 | % ------------------------------------------------------------------------------------------------------------- |
---|
920 | \subsection{Linear Bottom Friction (\protect\np{nn\_botfr} = 0 or 1) } |
---|
921 | \label{ZDF_bfr_linear} |
---|
922 | |
---|
923 | The linear bottom friction parameterisation (including the special case |
---|
924 | of a free-slip condition) assumes that the bottom friction |
---|
925 | is proportional to the interior velocity (i.e. the velocity of the last |
---|
926 | model level): |
---|
927 | \begin{equation} \label{Eq_zdfbfr_linear} |
---|
928 | {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \; \textbf{U}_h^b |
---|
929 | \end{equation} |
---|
930 | where $r$ is a friction coefficient expressed in ms$^{-1}$. |
---|
931 | This coefficient is generally estimated by setting a typical decay time |
---|
932 | $\tau$ in the deep ocean, |
---|
933 | and setting $r = H / \tau$, where $H$ is the ocean depth. Commonly accepted |
---|
934 | values of $\tau$ are of the order of 100 to 200 days \citep{Weatherly_JMR84}. |
---|
935 | A value $\tau^{-1} = 10^{-7}$~s$^{-1}$ equivalent to 115 days, is usually used |
---|
936 | in quasi-geostrophic models. One may consider the linear friction as an |
---|
937 | approximation of quadratic friction, $r \approx 2\;C_D\;U_{av}$ (\citet{Gill1982}, |
---|
938 | Eq. 9.6.6). For example, with a drag coefficient $C_D = 0.002$, a typical speed |
---|
939 | of tidal currents of $U_{av} =0.1$~m\;s$^{-1}$, and assuming an ocean depth |
---|
940 | $H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m\;s$^{-1}$. |
---|
941 | This is the default value used in \NEMO. It corresponds to a decay time scale |
---|
942 | of 115~days. It can be changed by specifying \np{rn\_bfri1} (namelist parameter). |
---|
943 | |
---|
944 | For the linear friction case the coefficients defined in the general |
---|
945 | expression \eqref{Eq_zdfbfr_bdef} are: |
---|
946 | \begin{equation} \label{Eq_zdfbfr_linbfr_b} |
---|
947 | \begin{split} |
---|
948 | c_b^u &= - r\\ |
---|
949 | c_b^v &= - r\\ |
---|
950 | \end{split} |
---|
951 | \end{equation} |
---|
952 | When \np{nn\_botfr}=1, the value of $r$ used is \np{rn\_bfri1}. |
---|
953 | Setting \np{nn\_botfr}=0 is equivalent to setting $r=0$ and leads to a free-slip |
---|
954 | bottom boundary condition. These values are assigned in \mdl{zdfbfr}. |
---|
955 | From v3.2 onwards there is support for local enhancement of these values |
---|
956 | via an externally defined 2D mask array (\np{ln\_bfr2d}=true) given |
---|
957 | in the \ifile{bfr\_coef} input NetCDF file. The mask values should vary from 0 to 1. |
---|
958 | Locations with a non-zero mask value will have the friction coefficient increased |
---|
959 | by $mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri1}. |
---|
960 | |
---|
961 | % ------------------------------------------------------------------------------------------------------------- |
---|
962 | % Non-Linear Bottom Friction |
---|
963 | % ------------------------------------------------------------------------------------------------------------- |
---|
964 | \subsection{Non-Linear Bottom Friction (\protect\np{nn\_botfr} = 2)} |
---|
965 | \label{ZDF_bfr_nonlinear} |
---|
966 | |
---|
967 | The non-linear bottom friction parameterisation assumes that the bottom |
---|
968 | friction is quadratic: |
---|
969 | \begin{equation} \label{Eq_zdfbfr_nonlinear} |
---|
970 | {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3 }\frac{\partial \textbf {U}_h |
---|
971 | }{\partial k}=C_D \;\sqrt {u_b ^2+v_b ^2+e_b } \;\; \textbf {U}_h^b |
---|
972 | \end{equation} |
---|
973 | where $C_D$ is a drag coefficient, and $e_b $ a bottom turbulent kinetic energy |
---|
974 | due to tides, internal waves breaking and other short time scale currents. |
---|
975 | A typical value of the drag coefficient is $C_D = 10^{-3} $. As an example, |
---|
976 | the CME experiment \citep{Treguier_JGR92} uses $C_D = 10^{-3}$ and |
---|
977 | $e_b = 2.5\;10^{-3}$m$^2$\;s$^{-2}$, while the FRAM experiment \citep{Killworth1992} |
---|
978 | uses $C_D = 1.4\;10^{-3}$ and $e_b =2.5\;\;10^{-3}$m$^2$\;s$^{-2}$. |
---|
979 | The CME choices have been set as default values (\np{rn\_bfri2} and \np{rn\_bfeb2} |
---|
980 | namelist parameters). |
---|
981 | |
---|
982 | As for the linear case, the bottom friction is imposed in the code by |
---|
983 | adding the trend due to the bottom friction to the general momentum trend |
---|
984 | in \mdl{dynbfr}. |
---|
985 | For the non-linear friction case the terms |
---|
986 | computed in \mdl{zdfbfr} are: |
---|
987 | \begin{equation} \label{Eq_zdfbfr_nonlinbfr} |
---|
988 | \begin{split} |
---|
989 | c_b^u &= - \; C_D\;\left[ u^2 + \left(\bar{\bar{v}}^{i+1,j}\right)^2 + e_b \right]^{1/2}\\ |
---|
990 | c_b^v &= - \; C_D\;\left[ \left(\bar{\bar{u}}^{i,j+1}\right)^2 + v^2 + e_b \right]^{1/2}\\ |
---|
991 | \end{split} |
---|
992 | \end{equation} |
---|
993 | |
---|
994 | The coefficients that control the strength of the non-linear bottom friction are |
---|
995 | initialised as namelist parameters: $C_D$= \np{rn\_bfri2}, and $e_b$ =\np{rn\_bfeb2}. |
---|
996 | Note for applications which treat tides explicitly a low or even zero value of |
---|
997 | \np{rn\_bfeb2} is recommended. From v3.2 onwards a local enhancement of $C_D$ is possible |
---|
998 | via an externally defined 2D mask array (\np{ln\_bfr2d}=true). This works in the same way |
---|
999 | as for the linear bottom friction case with non-zero masked locations increased by |
---|
1000 | $mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri2}. |
---|
1001 | |
---|
1002 | % ------------------------------------------------------------------------------------------------------------- |
---|
1003 | % Bottom Friction Log-layer |
---|
1004 | % ------------------------------------------------------------------------------------------------------------- |
---|
1005 | \subsection[Log-layer Bottom Friction enhancement (\protect\np{ln\_loglayer} = .true.)]{Log-layer Bottom Friction enhancement (\protect\np{nn\_botfr} = 2, \protect\np{ln\_loglayer} = .true.)} |
---|
1006 | \label{ZDF_bfr_loglayer} |
---|
1007 | |
---|
1008 | In the non-linear bottom friction case, the drag coefficient, $C_D$, can be optionally |
---|
1009 | enhanced using a "law of the wall" scaling. If \np{ln\_loglayer} = .true., $C_D$ is no |
---|
1010 | longer constant but is related to the thickness of the last wet layer in each column by: |
---|
1011 | |
---|
1012 | \begin{equation} |
---|
1013 | C_D = \left ( {\kappa \over {\rm log}\left ( 0.5e_{3t}/rn\_bfrz0 \right ) } \right )^2 |
---|
1014 | \end{equation} |
---|
1015 | |
---|
1016 | \noindent where $\kappa$ is the von-Karman constant and \np{rn\_bfrz0} is a roughness |
---|
1017 | length provided via the namelist. |
---|
1018 | |
---|
1019 | For stability, the drag coefficient is bounded such that it is kept greater or equal to |
---|
1020 | the base \np{rn\_bfri2} value and it is not allowed to exceed the value of an additional |
---|
1021 | namelist parameter: \np{rn\_bfri2\_max}, i.e.: |
---|
1022 | |
---|
1023 | \begin{equation} |
---|
1024 | rn\_bfri2 \leq C_D \leq rn\_bfri2\_max |
---|
1025 | \end{equation} |
---|
1026 | |
---|
1027 | \noindent Note also that a log-layer enhancement can also be applied to the top boundary |
---|
1028 | friction if under ice-shelf cavities are in use (\np{ln\_isfcav}=.true.). In this case, the |
---|
1029 | relevant namelist parameters are \np{rn\_tfrz0}, \np{rn\_tfri2} |
---|
1030 | and \np{rn\_tfri2\_max}. |
---|
1031 | |
---|
1032 | % ------------------------------------------------------------------------------------------------------------- |
---|
1033 | % Bottom Friction stability |
---|
1034 | % ------------------------------------------------------------------------------------------------------------- |
---|
1035 | \subsection{Bottom Friction stability considerations} |
---|
1036 | \label{ZDF_bfr_stability} |
---|
1037 | |
---|
1038 | Some care needs to exercised over the choice of parameters to ensure that the |
---|
1039 | implementation of bottom friction does not induce numerical instability. For |
---|
1040 | the purposes of stability analysis, an approximation to \eqref{Eq_zdfbfr_flux2} |
---|
1041 | is: |
---|
1042 | \begin{equation} \label{Eqn_bfrstab} |
---|
1043 | \begin{split} |
---|
1044 | \Delta u &= -\frac{{{\cal F}_h}^u}{e_{3u}}\;2 \rdt \\ |
---|
1045 | &= -\frac{ru}{e_{3u}}\;2\rdt\\ |
---|
1046 | \end{split} |
---|
1047 | \end{equation} |
---|
1048 | \noindent where linear bottom friction and a leapfrog timestep have been assumed. |
---|
1049 | To ensure that the bottom friction cannot reverse the direction of flow it is necessary to have: |
---|
1050 | \begin{equation} |
---|
1051 | |\Delta u| < \;|u| |
---|
1052 | \end{equation} |
---|
1053 | \noindent which, using \eqref{Eqn_bfrstab}, gives: |
---|
1054 | \begin{equation} |
---|
1055 | r\frac{2\rdt}{e_{3u}} < 1 \qquad \Rightarrow \qquad r < \frac{e_{3u}}{2\rdt}\\ |
---|
1056 | \end{equation} |
---|
1057 | This same inequality can also be derived in the non-linear bottom friction case |
---|
1058 | if a velocity of 1 m.s$^{-1}$ is assumed. Alternatively, this criterion can be |
---|
1059 | rearranged to suggest a minimum bottom box thickness to ensure stability: |
---|
1060 | \begin{equation} |
---|
1061 | e_{3u} > 2\;r\;\rdt |
---|
1062 | \end{equation} |
---|
1063 | \noindent which it may be necessary to impose if partial steps are being used. |
---|
1064 | For example, if $|u| = 1$ m.s$^{-1}$, $rdt = 1800$ s, $r = 10^{-3}$ then |
---|
1065 | $e_{3u}$ should be greater than 3.6 m. For most applications, with physically |
---|
1066 | sensible parameters these restrictions should not be of concern. But |
---|
1067 | caution may be necessary if attempts are made to locally enhance the bottom |
---|
1068 | friction parameters. |
---|
1069 | To ensure stability limits are imposed on the bottom friction coefficients both during |
---|
1070 | initialisation and at each time step. Checks at initialisation are made in \mdl{zdfbfr} |
---|
1071 | (assuming a 1 m.s$^{-1}$ velocity in the non-linear case). |
---|
1072 | The number of breaches of the stability criterion are reported as well as the minimum |
---|
1073 | and maximum values that have been set. The criterion is also checked at each time step, |
---|
1074 | using the actual velocity, in \mdl{dynbfr}. Values of the bottom friction coefficient are |
---|
1075 | reduced as necessary to ensure stability; these changes are not reported. |
---|
1076 | |
---|
1077 | Limits on the bottom friction coefficient are not imposed if the user has elected to |
---|
1078 | handle the bottom friction implicitly (see \S\ref{ZDF_bfr_imp}). The number of potential |
---|
1079 | breaches of the explicit stability criterion are still reported for information purposes. |
---|
1080 | |
---|
1081 | % ------------------------------------------------------------------------------------------------------------- |
---|
1082 | % Implicit Bottom Friction |
---|
1083 | % ------------------------------------------------------------------------------------------------------------- |
---|
1084 | \subsection[Implicit Bottom Friction (\protect\np{ln\_bfrimp})]{Implicit Bottom Friction (\protect\np{ln\_bfrimp}$=$\textit{T})} |
---|
1085 | \label{ZDF_bfr_imp} |
---|
1086 | |
---|
1087 | An optional implicit form of bottom friction has been implemented to improve |
---|
1088 | model stability. We recommend this option for shelf sea and coastal ocean applications, especially |
---|
1089 | for split-explicit time splitting. This option can be invoked by setting \np{ln\_bfrimp} |
---|
1090 | to \textit{true} in the \textit{nambfr} namelist. This option requires \np{ln\_zdfexp} to be \textit{false} |
---|
1091 | in the \textit{namzdf} namelist. |
---|
1092 | |
---|
1093 | This implementation is realised in \mdl{dynzdf\_imp} and \mdl{dynspg\_ts}. In \mdl{dynzdf\_imp}, the |
---|
1094 | bottom boundary condition is implemented implicitly. |
---|
1095 | |
---|
1096 | \begin{equation} \label{Eq_dynzdf_bfr} |
---|
1097 | \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{mbk} |
---|
1098 | = \binom{c_{b}^{u}u^{n+1}_{mbk}}{c_{b}^{v}v^{n+1}_{mbk}} |
---|
1099 | \end{equation} |
---|
1100 | |
---|
1101 | where $mbk$ is the layer number of the bottom wet layer. superscript $n+1$ means the velocity used in the |
---|
1102 | friction formula is to be calculated, so, it is implicit. |
---|
1103 | |
---|
1104 | If split-explicit time splitting is used, care must be taken to avoid the double counting of |
---|
1105 | the bottom friction in the 2-D barotropic momentum equations. As NEMO only updates the barotropic |
---|
1106 | pressure gradient and Coriolis' forcing terms in the 2-D barotropic calculation, we need to remove |
---|
1107 | the bottom friction induced by these two terms which has been included in the 3-D momentum trend |
---|
1108 | and update it with the latest value. On the other hand, the bottom friction contributed by the |
---|
1109 | other terms (e.g. the advection term, viscosity term) has been included in the 3-D momentum equations |
---|
1110 | and should not be added in the 2-D barotropic mode. |
---|
1111 | |
---|
1112 | The implementation of the implicit bottom friction in \mdl{dynspg\_ts} is done in two steps as the |
---|
1113 | following: |
---|
1114 | |
---|
1115 | \begin{equation} \label{Eq_dynspg_ts_bfr1} |
---|
1116 | \frac{\textbf{U}_{med}-\textbf{U}^{m-1}}{2\Delta t}=-g\nabla\eta-f\textbf{k}\times\textbf{U}^{m}+c_{b} |
---|
1117 | \left(\textbf{U}_{med}-\textbf{U}^{m-1}\right) |
---|
1118 | \end{equation} |
---|
1119 | \begin{equation} \label{Eq_dynspg_ts_bfr2} |
---|
1120 | \frac{\textbf{U}^{m+1}-\textbf{U}_{med}}{2\Delta t}=\textbf{T}+ |
---|
1121 | \left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{U}^{'}\right)- |
---|
1122 | 2\Delta t_{bc}c_{b}\left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{u}_{b}\right) |
---|
1123 | \end{equation} |
---|
1124 | |
---|
1125 | where $\textbf{T}$ is the vertical integrated 3-D momentum trend. We assume the leap-frog time-stepping |
---|
1126 | is used here. $\Delta t$ is the barotropic mode time step and $\Delta t_{bc}$ is the baroclinic mode time step. |
---|
1127 | $c_{b}$ is the friction coefficient. $\eta$ is the sea surface level calculated in the barotropic loops |
---|
1128 | while $\eta^{'}$ is the sea surface level used in the 3-D baroclinic mode. $\textbf{u}_{b}$ is the bottom |
---|
1129 | layer horizontal velocity. |
---|
1130 | |
---|
1131 | |
---|
1132 | |
---|
1133 | |
---|
1134 | % ------------------------------------------------------------------------------------------------------------- |
---|
1135 | % Bottom Friction with split-explicit time splitting |
---|
1136 | % ------------------------------------------------------------------------------------------------------------- |
---|
1137 | \subsection[Bottom Friction with split-explicit time splitting]{Bottom Friction with split-explicit time splitting (\protect\np{ln\_bfrimp})} |
---|
1138 | \label{ZDF_bfr_ts} |
---|
1139 | |
---|
1140 | When calculating the momentum trend due to bottom friction in \mdl{dynbfr}, the |
---|
1141 | bottom velocity at the before time step is used. This velocity includes both the |
---|
1142 | baroclinic and barotropic components which is appropriate when using either the |
---|
1143 | explicit or filtered surface pressure gradient algorithms (\key{dynspg\_exp} or |
---|
1144 | \key{dynspg\_flt}). Extra attention is required, however, when using |
---|
1145 | split-explicit time stepping (\key{dynspg\_ts}). In this case the free surface |
---|
1146 | equation is solved with a small time step \np{rn\_rdt}/\np{nn\_baro}, while the three |
---|
1147 | dimensional prognostic variables are solved with the longer time step |
---|
1148 | of \np{rn\_rdt} seconds. The trend in the barotropic momentum due to bottom |
---|
1149 | friction appropriate to this method is that given by the selected parameterisation |
---|
1150 | ($i.e.$ linear or non-linear bottom friction) computed with the evolving velocities |
---|
1151 | at each barotropic timestep. |
---|
1152 | |
---|
1153 | In the case of non-linear bottom friction, we have elected to partially linearise |
---|
1154 | the problem by keeping the coefficients fixed throughout the barotropic |
---|
1155 | time-stepping to those computed in \mdl{zdfbfr} using the now timestep. |
---|
1156 | This decision allows an efficient use of the $c_b^{\vect{U}}$ coefficients to: |
---|
1157 | |
---|
1158 | \begin{enumerate} |
---|
1159 | \item On entry to \rou{dyn\_spg\_ts}, remove the contribution of the before |
---|
1160 | barotropic velocity to the bottom friction component of the vertically |
---|
1161 | integrated momentum trend. Note the same stability check that is carried out |
---|
1162 | on the bottom friction coefficient in \mdl{dynbfr} has to be applied here to |
---|
1163 | ensure that the trend removed matches that which was added in \mdl{dynbfr}. |
---|
1164 | \item At each barotropic step, compute the contribution of the current barotropic |
---|
1165 | velocity to the trend due to bottom friction. Add this contribution to the |
---|
1166 | vertically integrated momentum trend. This contribution is handled implicitly which |
---|
1167 | eliminates the need to impose a stability criteria on the values of the bottom friction |
---|
1168 | coefficient within the barotropic loop. |
---|
1169 | \end{enumerate} |
---|
1170 | |
---|
1171 | Note that the use of an implicit formulation within the barotropic loop |
---|
1172 | for the bottom friction trend means that any limiting of the bottom friction coefficient |
---|
1173 | in \mdl{dynbfr} does not adversely affect the solution when using split-explicit time |
---|
1174 | splitting. This is because the major contribution to bottom friction is likely to come from |
---|
1175 | the barotropic component which uses the unrestricted value of the coefficient. However, if the |
---|
1176 | limiting is thought to be having a major effect (a more likely prospect in coastal and shelf seas |
---|
1177 | applications) then the fully implicit form of the bottom friction should be used (see \S\ref{ZDF_bfr_imp} ) |
---|
1178 | which can be selected by setting \np{ln\_bfrimp} $=$ \textit{true}. |
---|
1179 | |
---|
1180 | Otherwise, the implicit formulation takes the form: |
---|
1181 | \begin{equation} \label{Eq_zdfbfr_implicitts} |
---|
1182 | \bar{U}^{t+ \rdt} = \; \left [ \bar{U}^{t-\rdt}\; + 2 \rdt\;RHS \right ] / \left [ 1 - 2 \rdt \;c_b^{u} / H_e \right ] |
---|
1183 | \end{equation} |
---|
1184 | where $\bar U$ is the barotropic velocity, $H_e$ is the full depth (including sea surface height), |
---|
1185 | $c_b^u$ is the bottom friction coefficient as calculated in \rou{zdf\_bfr} and $RHS$ represents |
---|
1186 | all the components to the vertically integrated momentum trend except for that due to bottom friction. |
---|
1187 | |
---|
1188 | |
---|
1189 | |
---|
1190 | |
---|
1191 | % ================================================================ |
---|
1192 | % Tidal Mixing |
---|
1193 | % ================================================================ |
---|
1194 | \section{Tidal Mixing (\protect\key{zdftmx})} |
---|
1195 | \label{ZDF_tmx} |
---|
1196 | |
---|
1197 | %--------------------------------------------namzdf_tmx-------------------------------------------------- |
---|
1198 | \namdisplay{namzdf_tmx} |
---|
1199 | %-------------------------------------------------------------------------------------------------------------- |
---|
1200 | |
---|
1201 | |
---|
1202 | % ------------------------------------------------------------------------------------------------------------- |
---|
1203 | % Bottom intensified tidal mixing |
---|
1204 | % ------------------------------------------------------------------------------------------------------------- |
---|
1205 | \subsection{Bottom intensified tidal mixing} |
---|
1206 | \label{ZDF_tmx_bottom} |
---|
1207 | |
---|
1208 | Options are defined through the \ngn{namzdf\_tmx} namelist variables. |
---|
1209 | The parameterization of tidal mixing follows the general formulation for |
---|
1210 | the vertical eddy diffusivity proposed by \citet{St_Laurent_al_GRL02} and |
---|
1211 | first introduced in an OGCM by \citep{Simmons_al_OM04}. |
---|
1212 | In this formulation an additional vertical diffusivity resulting from internal tide breaking, |
---|
1213 | $A^{vT}_{tides}$ is expressed as a function of $E(x,y)$, the energy transfer from barotropic |
---|
1214 | tides to baroclinic tides : |
---|
1215 | \begin{equation} \label{Eq_Ktides} |
---|
1216 | A^{vT}_{tides} = q \,\Gamma \,\frac{ E(x,y) \, F(z) }{ \rho \, N^2 } |
---|
1217 | \end{equation} |
---|
1218 | where $\Gamma$ is the mixing efficiency, $N$ the Brunt-Vais\"{a}l\"{a} frequency |
---|
1219 | (see \S\ref{TRA_bn2}), $\rho$ the density, $q$ the tidal dissipation efficiency, |
---|
1220 | and $F(z)$ the vertical structure function. |
---|
1221 | |
---|
1222 | The mixing efficiency of turbulence is set by $\Gamma$ (\np{rn\_me} namelist parameter) |
---|
1223 | and is usually taken to be the canonical value of $\Gamma = 0.2$ (Osborn 1980). |
---|
1224 | The tidal dissipation efficiency is given by the parameter $q$ (\np{rn\_tfe} namelist parameter) |
---|
1225 | represents the part of the internal wave energy flux $E(x, y)$ that is dissipated locally, |
---|
1226 | with the remaining $1-q$ radiating away as low mode internal waves and |
---|
1227 | contributing to the background internal wave field. A value of $q=1/3$ is typically used |
---|
1228 | \citet{St_Laurent_al_GRL02}. |
---|
1229 | The vertical structure function $F(z)$ models the distribution of the turbulent mixing in the vertical. |
---|
1230 | It is implemented as a simple exponential decaying upward away from the bottom, |
---|
1231 | with a vertical scale of $h_o$ (\np{rn\_htmx} namelist parameter, with a typical value of $500\,m$) \citep{St_Laurent_Nash_DSR04}, |
---|
1232 | \begin{equation} \label{Eq_Fz} |
---|
1233 | F(i,j,k) = \frac{ e^{ -\frac{H+z}{h_o} } }{ h_o \left( 1- e^{ -\frac{H}{h_o} } \right) } |
---|
1234 | \end{equation} |
---|
1235 | and is normalized so that vertical integral over the water column is unity. |
---|
1236 | |
---|
1237 | The associated vertical viscosity is calculated from the vertical |
---|
1238 | diffusivity assuming a Prandtl number of 1, $i.e.$ $A^{vm}_{tides}=A^{vT}_{tides}$. |
---|
1239 | In the limit of $N \rightarrow 0$ (or becoming negative), the vertical diffusivity |
---|
1240 | is capped at $300\,cm^2/s$ and impose a lower limit on $N^2$ of \np{rn\_n2min} |
---|
1241 | usually set to $10^{-8} s^{-2}$. These bounds are usually rarely encountered. |
---|
1242 | |
---|
1243 | The internal wave energy map, $E(x, y)$ in \eqref{Eq_Ktides}, is derived |
---|
1244 | from a barotropic model of the tides utilizing a parameterization of the |
---|
1245 | conversion of barotropic tidal energy into internal waves. |
---|
1246 | The essential goal of the parameterization is to represent the momentum |
---|
1247 | exchange between the barotropic tides and the unrepresented internal waves |
---|
1248 | induced by the tidal flow over rough topography in a stratified ocean. |
---|
1249 | In the current version of \NEMO, the map is built from the output of |
---|
1250 | the barotropic global ocean tide model MOG2D-G \citep{Carrere_Lyard_GRL03}. |
---|
1251 | This model provides the dissipation associated with internal wave energy for the M2 and K1 |
---|
1252 | tides component (Fig.~\ref{Fig_ZDF_M2_K1_tmx}). The S2 dissipation is simply approximated |
---|
1253 | as being $1/4$ of the M2 one. The internal wave energy is thus : $E(x, y) = 1.25 E_{M2} + E_{K1}$. |
---|
1254 | Its global mean value is $1.1$ TW, in agreement with independent estimates |
---|
1255 | \citep{Egbert_Ray_Nat00, Egbert_Ray_JGR01}. |
---|
1256 | |
---|
1257 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
1258 | \begin{figure}[!t] \begin{center} |
---|
1259 | \includegraphics[width=0.90\textwidth]{Fig_ZDF_M2_K1_tmx} |
---|
1260 | \caption{ \protect\label{Fig_ZDF_M2_K1_tmx} |
---|
1261 | (a) M2 and (b) K1 internal wave drag energy from \citet{Carrere_Lyard_GRL03} ($W/m^2$). } |
---|
1262 | \end{center} \end{figure} |
---|
1263 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
1264 | |
---|
1265 | % ------------------------------------------------------------------------------------------------------------- |
---|
1266 | % Indonesian area specific treatment |
---|
1267 | % ------------------------------------------------------------------------------------------------------------- |
---|
1268 | \subsection{Indonesian area specific treatment (\protect\np{ln\_zdftmx\_itf})} |
---|
1269 | \label{ZDF_tmx_itf} |
---|
1270 | |
---|
1271 | When the Indonesian Through Flow (ITF) area is included in the model domain, |
---|
1272 | a specific treatment of tidal induced mixing in this area can be used. |
---|
1273 | It is activated through the namelist logical \np{ln\_tmx\_itf}, and the user must provide |
---|
1274 | an input NetCDF file, \ifile{mask\_itf}, which contains a mask array defining the ITF area |
---|
1275 | where the specific treatment is applied. |
---|
1276 | |
---|
1277 | When \np{ln\_tmx\_itf}=true, the two key parameters $q$ and $F(z)$ are adjusted following |
---|
1278 | the parameterisation developed by \citet{Koch-Larrouy_al_GRL07}: |
---|
1279 | |
---|
1280 | First, the Indonesian archipelago is a complex geographic region |
---|
1281 | with a series of large, deep, semi-enclosed basins connected via |
---|
1282 | numerous narrow straits. Once generated, internal tides remain |
---|
1283 | confined within this semi-enclosed area and hardly radiate away. |
---|
1284 | Therefore all the internal tides energy is consumed within this area. |
---|
1285 | So it is assumed that $q = 1$, $i.e.$ all the energy generated is available for mixing. |
---|
1286 | Note that for test purposed, the ITF tidal dissipation efficiency is a |
---|
1287 | namelist parameter (\np{rn\_tfe\_itf}). A value of $1$ or close to is |
---|
1288 | this recommended for this parameter. |
---|
1289 | |
---|
1290 | Second, the vertical structure function, $F(z)$, is no more associated |
---|
1291 | with a bottom intensification of the mixing, but with a maximum of |
---|
1292 | energy available within the thermocline. \citet{Koch-Larrouy_al_GRL07} |
---|
1293 | have suggested that the vertical distribution of the energy dissipation |
---|
1294 | proportional to $N^2$ below the core of the thermocline and to $N$ above. |
---|
1295 | The resulting $F(z)$ is: |
---|
1296 | \begin{equation} \label{Eq_Fz_itf} |
---|
1297 | F(i,j,k) \sim \left\{ \begin{aligned} |
---|
1298 | \frac{q\,\Gamma E(i,j) } {\rho N \, \int N dz} \qquad \text{when $\partial_z N < 0$} \\ |
---|
1299 | \frac{q\,\Gamma E(i,j) } {\rho \, \int N^2 dz} \qquad \text{when $\partial_z N > 0$} |
---|
1300 | \end{aligned} \right. |
---|
1301 | \end{equation} |
---|
1302 | |
---|
1303 | Averaged over the ITF area, the resulting tidal mixing coefficient is $1.5\,cm^2/s$, |
---|
1304 | which agrees with the independent estimates inferred from observations. |
---|
1305 | Introduced in a regional OGCM, the parameterization improves the water mass |
---|
1306 | characteristics in the different Indonesian seas, suggesting that the horizontal |
---|
1307 | and vertical distributions of the mixing are adequately prescribed |
---|
1308 | \citep{Koch-Larrouy_al_GRL07, Koch-Larrouy_al_OD08a, Koch-Larrouy_al_OD08b}. |
---|
1309 | Note also that such a parameterisation has a significant impact on the behaviour |
---|
1310 | of global coupled GCMs \citep{Koch-Larrouy_al_CD10}. |
---|
1311 | |
---|
1312 | |
---|
1313 | % ================================================================ |
---|
1314 | % Internal wave-driven mixing |
---|
1315 | % ================================================================ |
---|
1316 | \section{Internal wave-driven mixing (\protect\key{zdftmx\_new})} |
---|
1317 | \label{ZDF_tmx_new} |
---|
1318 | |
---|
1319 | %--------------------------------------------namzdf_tmx_new------------------------------------------ |
---|
1320 | \namdisplay{namzdf_tmx_new} |
---|
1321 | %-------------------------------------------------------------------------------------------------------------- |
---|
1322 | |
---|
1323 | The parameterization of mixing induced by breaking internal waves is a generalization |
---|
1324 | of the approach originally proposed by \citet{St_Laurent_al_GRL02}. |
---|
1325 | A three-dimensional field of internal wave energy dissipation $\epsilon(x,y,z)$ is first constructed, |
---|
1326 | and the resulting diffusivity is obtained as |
---|
1327 | \begin{equation} \label{Eq_Kwave} |
---|
1328 | A^{vT}_{wave} = R_f \,\frac{ \epsilon }{ \rho \, N^2 } |
---|
1329 | \end{equation} |
---|
1330 | where $R_f$ is the mixing efficiency and $\epsilon$ is a specified three dimensional distribution |
---|
1331 | of the energy available for mixing. If the \np{ln\_mevar} namelist parameter is set to false, |
---|
1332 | the mixing efficiency is taken as constant and equal to 1/6 \citep{Osborn_JPO80}. |
---|
1333 | In the opposite (recommended) case, $R_f$ is instead a function of the turbulence intensity parameter |
---|
1334 | $Re_b = \frac{ \epsilon}{\nu \, N^2}$, with $\nu$ the molecular viscosity of seawater, |
---|
1335 | following the model of \cite{Bouffard_Boegman_DAO2013} |
---|
1336 | and the implementation of \cite{de_lavergne_JPO2016_efficiency}. |
---|
1337 | Note that $A^{vT}_{wave}$ is bounded by $10^{-2}\,m^2/s$, a limit that is often reached when the mixing efficiency is constant. |
---|
1338 | |
---|
1339 | In addition to the mixing efficiency, the ratio of salt to heat diffusivities can chosen to vary |
---|
1340 | as a function of $Re_b$ by setting the \np{ln\_tsdiff} parameter to true, a recommended choice). |
---|
1341 | This parameterization of differential mixing, due to \cite{Jackson_Rehmann_JPO2014}, |
---|
1342 | is implemented as in \cite{de_lavergne_JPO2016_efficiency}. |
---|
1343 | |
---|
1344 | The three-dimensional distribution of the energy available for mixing, $\epsilon(i,j,k)$, is constructed |
---|
1345 | from three static maps of column-integrated internal wave energy dissipation, $E_{cri}(i,j)$, |
---|
1346 | $E_{pyc}(i,j)$, and $E_{bot}(i,j)$, combined to three corresponding vertical structures |
---|
1347 | (de Lavergne et al., in prep): |
---|
1348 | \begin{align*} |
---|
1349 | F_{cri}(i,j,k) &\propto e^{-h_{ab} / h_{cri} }\\ |
---|
1350 | F_{pyc}(i,j,k) &\propto N^{n\_p}\\ |
---|
1351 | F_{bot}(i,j,k) &\propto N^2 \, e^{- h_{wkb} / h_{bot} } |
---|
1352 | \end{align*} |
---|
1353 | In the above formula, $h_{ab}$ denotes the height above bottom, |
---|
1354 | $h_{wkb}$ denotes the WKB-stretched height above bottom, defined by |
---|
1355 | \begin{equation*} |
---|
1356 | h_{wkb} = H \, \frac{ \int_{-H}^{z} N \, dz' } { \int_{-H}^{\eta} N \, dz' } \; , |
---|
1357 | \end{equation*} |
---|
1358 | The $n_p$ parameter (given by \np{nn\_zpyc} in \ngn{namzdf\_tmx\_new} namelist) controls the stratification-dependence of the pycnocline-intensified dissipation. |
---|
1359 | It can take values of 1 (recommended) or 2. |
---|
1360 | Finally, the vertical structures $F_{cri}$ and $F_{bot}$ require the specification of |
---|
1361 | the decay scales $h_{cri}(i,j)$ and $h_{bot}(i,j)$, which are defined by two additional input maps. |
---|
1362 | $h_{cri}$ is related to the large-scale topography of the ocean (etopo2) |
---|
1363 | and $h_{bot}$ is a function of the energy flux $E_{bot}$, the characteristic horizontal scale of |
---|
1364 | the abyssal hill topography \citep{Goff_JGR2010} and the latitude. |
---|
1365 | |
---|
1366 | % ================================================================ |
---|
1367 | |
---|
1368 | |
---|
1369 | |
---|
1370 | \end{document} |
---|