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