[10414] | 1 | \documentclass[../main/NEMO_manual]{subfiles} |
---|
| 2 | |
---|
[6997] | 3 | \begin{document} |
---|
[10501] | 4 | |
---|
[9393] | 5 | \chapter{Model Basics} |
---|
[13463] | 6 | \label{chap:MB} |
---|
[707] | 7 | |
---|
[13463] | 8 | \thispagestyle{plain} |
---|
[2282] | 9 | |
---|
[13463] | 10 | \chaptertoc |
---|
| 11 | |
---|
| 12 | \paragraph{Changes record} ~\\ |
---|
| 13 | |
---|
| 14 | {\footnotesize |
---|
| 15 | \begin{tabular}{l||l|l} |
---|
| 16 | Release & Author(s) & Modifications \\ |
---|
| 17 | \hline |
---|
| 18 | {\em 4.0} & {\em Mike Bell } & {\em Review } \\ |
---|
| 19 | {\em 3.6} & {\em Tim Graham and Gurvan Madec } & {\em Updates } \\ |
---|
| 20 | {\em $\leq$ 3.4} & {\em Gurvan Madec and S\'{e}bastien Masson} & {\em First version} \\ |
---|
| 21 | \end{tabular} |
---|
| 22 | } |
---|
| 23 | |
---|
| 24 | \clearpage |
---|
| 25 | |
---|
| 26 | %% ================================================================================================= |
---|
[9393] | 27 | \section{Primitive equations} |
---|
[13463] | 28 | \label{sec:MB_PE} |
---|
[707] | 29 | |
---|
[13463] | 30 | %% ================================================================================================= |
---|
[9393] | 31 | \subsection{Vector invariant formulation} |
---|
[13463] | 32 | \label{subsec:MB_PE_vector} |
---|
[707] | 33 | |
---|
[10354] | 34 | The ocean is a fluid that can be described to a good approximation by the primitive equations, |
---|
[13463] | 35 | \ie\ the Navier-Stokes equations along with a nonlinear equation of state which |
---|
[10354] | 36 | couples the two active tracers (temperature and salinity) to the fluid velocity, |
---|
| 37 | plus the following additional assumptions made from scale considerations: |
---|
[707] | 38 | |
---|
[13463] | 39 | \begin{labeling}{Neglect of additional Coriolis terms} |
---|
| 40 | \item [\textit{Spherical Earth approximation}] The geopotential surfaces are assumed to |
---|
| 41 | be oblate spheroids that follow the Earth's bulge; |
---|
| 42 | these spheroids are approximated by spheres with gravity locally vertical |
---|
| 43 | (parallel to the Earth's radius) and independent of latitude |
---|
| 44 | \citep[][section 2]{white.hoskins.ea_QJRMS05}. |
---|
| 45 | \item [\textit{Thin-shell approximation}] The ocean depth is neglected compared to the earth's radius |
---|
| 46 | \item [\textit{Turbulent closure hypothesis}] The turbulent fluxes |
---|
[10501] | 47 | (which represent the effect of small scale processes on the large-scale) |
---|
| 48 | are expressed in terms of large-scale features |
---|
[13463] | 49 | \item [\textit{Boussinesq hypothesis}] Density variations are neglected except in |
---|
| 50 | their contribution to the buoyancy force |
---|
[10414] | 51 | \begin{equation} |
---|
[13463] | 52 | \label{eq:MB_PE_eos} |
---|
[10501] | 53 | \rho = \rho \ (T,S,p) |
---|
[707] | 54 | \end{equation} |
---|
[13463] | 55 | \item [\textit{Hydrostatic hypothesis}] The vertical momentum equation is reduced to |
---|
| 56 | a balance between the vertical pressure gradient and the buoyancy force |
---|
[10501] | 57 | (this removes convective processes from the initial Navier-Stokes equations and so |
---|
| 58 | convective processes must be parameterized instead) |
---|
[10414] | 59 | \begin{equation} |
---|
[13463] | 60 | \label{eq:MB_PE_hydrostatic} |
---|
[10501] | 61 | \pd[p]{z} = - \rho \ g |
---|
[707] | 62 | \end{equation} |
---|
[13463] | 63 | \item [\textit{Incompressibility hypothesis}] The three dimensional divergence of |
---|
| 64 | the velocity vector $\vect U$ is assumed to be zero. |
---|
[10414] | 65 | \begin{equation} |
---|
[13463] | 66 | \label{eq:MB_PE_continuity} |
---|
[10501] | 67 | \nabla \cdot \vect U = 0 |
---|
[707] | 68 | \end{equation} |
---|
[13463] | 69 | \item [\textit{Neglect of additional Coriolis terms}] The Coriolis terms that vary with |
---|
| 70 | the cosine of latitude are neglected. |
---|
| 71 | These terms may be non-negligible where the Brunt-V\"{a}is\"{a}l\"{a} frequency $N$ is small, |
---|
| 72 | either in the deep ocean or in the sub-mesoscale motions of the mixed layer, |
---|
| 73 | or near the equator \citep[][section 1]{white.hoskins.ea_QJRMS05}. |
---|
| 74 | They can be consistently included as part of the ocean dynamics |
---|
| 75 | \citep[][section 3(d)]{white.hoskins.ea_QJRMS05} and are retained in the MIT ocean model. |
---|
| 76 | \end{labeling} |
---|
[10501] | 77 | |
---|
| 78 | Because the gravitational force is so dominant in the equations of large-scale motions, |
---|
[11335] | 79 | it is useful to choose an orthogonal set of unit vectors $(i,j,k)$ linked to the Earth such that |
---|
[10501] | 80 | $k$ is the local upward vector and $(i,j)$ are two vectors orthogonal to $k$, |
---|
[13463] | 81 | \ie\ tangent to the geopotential surfaces. |
---|
| 82 | Let us define the following variables: |
---|
| 83 | $\vect U$ the vector velocity, $\vect U = \vect U_h + w \, \vect k$ |
---|
| 84 | (the subscript $h$ denotes the local horizontal vector, \ie\ over the $(i,j)$ plane), |
---|
[10501] | 85 | $T$ the potential temperature, $S$ the salinity, $\rho$ the \textit{in situ} density. |
---|
| 86 | The vector invariant form of the primitive equations in the $(i,j,k)$ vector system provides |
---|
| 87 | the following equations: |
---|
| 88 | \begin{subequations} |
---|
[13463] | 89 | \label{eq:MB_PE} |
---|
[10501] | 90 | \begin{gather} |
---|
[13463] | 91 | \shortintertext{$-$ the momentum balance} |
---|
| 92 | \label{eq:MB_PE_dyn} |
---|
| 93 | \pd[\vect U_h]{t} = - \lt[ (\nabla \times \vect U) \times \vect U + \frac{1}{2} \nabla \lt( \vect U^2 \rt) \rt]_h - f \; k \times \vect U_h - \frac{1}{\rho_o} \nabla_h p + \vect D^{\vect U} + \vect F^{\vect U} |
---|
| 94 | \shortintertext{$-$ the heat and salt conservation equations} |
---|
| 95 | \label{eq:MB_PE_tra_T} |
---|
[10501] | 96 | \pd[T]{t} = - \nabla \cdot (T \ \vect U) + D^T + F^T \\ |
---|
[13463] | 97 | \label{eq:MB_PE_tra_S} |
---|
[10501] | 98 | \pd[S]{t} = - \nabla \cdot (S \ \vect U) + D^S + F^S |
---|
| 99 | \end{gather} |
---|
[707] | 100 | \end{subequations} |
---|
[10501] | 101 | where $\nabla$ is the generalised derivative vector operator in $(i,j,k)$ directions, $t$ is the time, |
---|
| 102 | $z$ is the vertical coordinate, $\rho$ is the \textit{in situ} density given by the equation of state |
---|
[13463] | 103 | (\autoref{eq:MB_PE_eos}), $\rho_o$ is a reference density, $p$ the pressure, |
---|
[10501] | 104 | $f = 2 \vect \Omega \cdot k$ is the Coriolis acceleration |
---|
[13463] | 105 | (where $\vect \Omega$ is the Earth's angular velocity vector), |
---|
| 106 | and $g$ is the gravitational acceleration. |
---|
[10501] | 107 | $\vect D^{\vect U}$, $D^T$ and $D^S$ are the parameterisations of small-scale physics for momentum, |
---|
| 108 | temperature and salinity, and $\vect F^{\vect U}$, $F^T$ and $F^S$ surface forcing terms. |
---|
[13463] | 109 | Their nature and formulation are discussed in \autoref{sec:MB_zdf_ldf} and |
---|
| 110 | \autoref{subsec:MB_boundary_condition}. |
---|
[707] | 111 | |
---|
[13463] | 112 | %% ================================================================================================= |
---|
[9393] | 113 | \subsection{Boundary conditions} |
---|
[13463] | 114 | \label{subsec:MB_boundary_condition} |
---|
[707] | 115 | |
---|
[10354] | 116 | An ocean is bounded by complex coastlines, bottom topography at its base and |
---|
| 117 | an air-sea or ice-sea interface at its top. |
---|
[10501] | 118 | These boundaries can be defined by two surfaces, $z = - H(i,j)$ and $z = \eta(i,j,k,t)$, |
---|
[13463] | 119 | where $H$ is the depth of the ocean bottom and $\eta$ is the height of the sea surface |
---|
| 120 | (discretisation can introduce additional artificial ``side-wall'' boundaries). |
---|
| 121 | Both $H$ and $\eta$ are referenced to a surface of constant geopotential |
---|
| 122 | (\ie\ a mean sea surface height) on which $z = 0$ (\autoref{fig:MB_ocean_bc}). |
---|
| 123 | Through these two boundaries, the ocean can exchange fluxes of heat, fresh water, salt, |
---|
| 124 | and momentum with the solid earth, the continental margins, the sea ice and the atmosphere. |
---|
| 125 | However, some of these fluxes are so weak that |
---|
| 126 | even on climatic time scales of thousands of years they can be neglected. |
---|
[10354] | 127 | In the following, we briefly review the fluxes exchanged at the interfaces between the ocean and |
---|
| 128 | the other components of the earth system. |
---|
[707] | 129 | |
---|
[13463] | 130 | \begin{figure} |
---|
| 131 | \centering |
---|
| 132 | \includegraphics[width=0.66\textwidth]{MB_ocean_bc} |
---|
| 133 | \caption[Ocean boundary conditions]{ |
---|
| 134 | The ocean is bounded by two surfaces, $z = - H(i,j)$ and $z = \eta(i,j,t)$, |
---|
| 135 | where $H$ is the depth of the sea floor and $\eta$ the height of the sea surface. |
---|
| 136 | Both $H$ and $\eta$ are referenced to $z = 0$.} |
---|
| 137 | \label{fig:MB_ocean_bc} |
---|
[10354] | 138 | \end{figure} |
---|
[707] | 139 | |
---|
| 140 | \begin{description} |
---|
[13463] | 141 | \item [Land - ocean] The major flux between continental margins and the ocean is a mass exchange of |
---|
| 142 | fresh water through river runoff. |
---|
[10354] | 143 | Such an exchange modifies the sea surface salinity especially in the vicinity of major river mouths. |
---|
[13463] | 144 | It can be neglected for short range integrations but |
---|
| 145 | has to be taken into account for long term integrations as |
---|
[10354] | 146 | it influences the characteristics of water masses formed (especially at high latitudes). |
---|
| 147 | It is required in order to close the water cycle of the climate system. |
---|
[13463] | 148 | It is usually specified as a fresh water flux at the air-sea interface in |
---|
| 149 | the vicinity of river mouths. |
---|
| 150 | \item [Solid earth - ocean] Heat and salt fluxes through the sea floor are small, |
---|
| 151 | except in special areas of little extent. |
---|
| 152 | They are usually neglected in the model \footnote{ |
---|
[10354] | 153 | In fact, it has been shown that the heat flux associated with the solid Earth cooling |
---|
[13463] | 154 | (\ie\ the geothermal heating) is not negligible for |
---|
| 155 | the thermohaline circulation of the world ocean (see \autoref{subsec:TRA_bbc}). |
---|
[10354] | 156 | }. |
---|
| 157 | The boundary condition is thus set to no flux of heat and salt across solid boundaries. |
---|
[13463] | 158 | For momentum, the situation is different. |
---|
| 159 | There is no flow across solid boundaries, |
---|
| 160 | \ie\ the velocity normal to the ocean bottom and coastlines is zero |
---|
| 161 | (in other words, the bottom velocity is parallel to solid boundaries). |
---|
| 162 | This kinematic boundary condition can be expressed as: |
---|
[10414] | 163 | \begin{equation} |
---|
[13463] | 164 | \label{eq:MB_w_bbc} |
---|
[10501] | 165 | w = - \vect U_h \cdot \nabla_h (H) |
---|
[10354] | 166 | \end{equation} |
---|
| 167 | In addition, the ocean exchanges momentum with the earth through frictional processes. |
---|
| 168 | Such momentum transfer occurs at small scales in a boundary layer. |
---|
[13463] | 169 | It must be parameterized in terms of turbulent fluxes using |
---|
| 170 | bottom and/or lateral boundary conditions. |
---|
[10354] | 171 | Its specification depends on the nature of the physical parameterisation used for |
---|
[13463] | 172 | $\vect D^{\vect U}$ in \autoref{eq:MB_PE_dyn}. |
---|
| 173 | It is discussed in \autoref{eq:MB_zdf}. % and Chap. III.6 to 9. |
---|
| 174 | \item [Atmosphere - ocean] The kinematic surface condition plus the mass flux of fresh water PE |
---|
| 175 | (the precipitation minus evaporation budget) leads to: |
---|
[10414] | 176 | \[ |
---|
[13463] | 177 | % \label{eq:MB_w_sbc} |
---|
[10501] | 178 | w = \pd[\eta]{t} + \lt. \vect U_h \rt|_{z = \eta} \cdot \nabla_h (\eta) + P - E |
---|
[10414] | 179 | \] |
---|
[13463] | 180 | The dynamic boundary condition, neglecting the surface tension |
---|
| 181 | (which removes capillary waves from the system) leads to |
---|
| 182 | the continuity of pressure across the interface $z = \eta$. |
---|
[10354] | 183 | The atmosphere and ocean also exchange horizontal momentum (wind stress), and heat. |
---|
[13463] | 184 | \item [Sea ice - ocean] The ocean and sea ice exchange heat, salt, fresh water and momentum. |
---|
[10354] | 185 | The sea surface temperature is constrained to be at the freezing point at the interface. |
---|
[10501] | 186 | Sea ice salinity is very low ($\sim4-6 \, psu$) compared to those of the ocean ($\sim34 \, psu$). |
---|
[13463] | 187 | The cycle of freezing/melting is associated with fresh water and salt fluxes that |
---|
| 188 | cannot be neglected. |
---|
[707] | 189 | \end{description} |
---|
| 190 | |
---|
[13463] | 191 | %% ================================================================================================= |
---|
[10501] | 192 | \section{Horizontal pressure gradient} |
---|
[13463] | 193 | \label{sec:MB_hor_pg} |
---|
[707] | 194 | |
---|
[13463] | 195 | %% ================================================================================================= |
---|
[9393] | 196 | \subsection{Pressure formulation} |
---|
[13463] | 197 | \label{subsec:MB_p_formulation} |
---|
[707] | 198 | |
---|
[10354] | 199 | The total pressure at a given depth $z$ is composed of a surface pressure $p_s$ at |
---|
[10501] | 200 | a reference geopotential surface ($z = 0$) and a hydrostatic pressure $p_h$ such that: |
---|
| 201 | $p(i,j,k,t) = p_s(i,j,t) + p_h(i,j,k,t)$. |
---|
[13463] | 202 | The latter is computed by integrating (\autoref{eq:MB_PE_hydrostatic}), |
---|
| 203 | assuming that pressure in decibars can be approximated by depth in meters in (\autoref{eq:MB_PE_eos}). |
---|
[1224] | 204 | The hydrostatic pressure is then given by: |
---|
[10414] | 205 | \[ |
---|
[13463] | 206 | % \label{eq:MB_pressure} |
---|
[10501] | 207 | p_h (i,j,z,t) = \int_{\varsigma = z}^{\varsigma = 0} g \; \rho (T,S,\varsigma) \; d \varsigma |
---|
[10414] | 208 | \] |
---|
[10354] | 209 | Two strategies can be considered for the surface pressure term: |
---|
[13463] | 210 | \begin{enumerate*}[label=(\textit{\alph*})] |
---|
| 211 | \item introduce of a new variable $\eta$, the free-surface elevation, |
---|
[10354] | 212 | for which a prognostic equation can be established and solved; |
---|
[13463] | 213 | \item assume that the ocean surface is a rigid lid, |
---|
[10354] | 214 | on which the pressure (or its horizontal gradient) can be diagnosed. |
---|
[13463] | 215 | \end{enumerate*} |
---|
[10354] | 216 | When the former strategy is used, one solution of the free-surface elevation consists of |
---|
| 217 | the excitation of external gravity waves. |
---|
| 218 | The flow is barotropic and the surface moves up and down with gravity as the restoring force. |
---|
| 219 | The phase speed of such waves is high (some hundreds of metres per second) so that |
---|
[11335] | 220 | the time step has to be very short when they are present in the model. |
---|
[10501] | 221 | The latter strategy filters out these waves since the rigid lid approximation implies $\eta = 0$, |
---|
[13463] | 222 | \ie\ the sea surface is the surface $z = 0$. |
---|
[10354] | 223 | This well known approximation increases the surface wave speed to infinity and |
---|
[13463] | 224 | modifies certain other longwave dynamics (\eg\ barotropic Rossby or planetary waves). |
---|
[10354] | 225 | The rigid-lid hypothesis is an obsolescent feature in modern OGCMs. |
---|
[13463] | 226 | It has been available until the release 3.1 of \NEMO, |
---|
| 227 | and it has been removed in release 3.2 and followings. |
---|
[11335] | 228 | Only the free surface formulation is now described in this document (see the next sub-section). |
---|
[707] | 229 | |
---|
[13463] | 230 | %% ================================================================================================= |
---|
[9393] | 231 | \subsection{Free surface formulation} |
---|
[13463] | 232 | \label{subsec:MB_free_surface} |
---|
[707] | 233 | |
---|
[10354] | 234 | In the free surface formulation, a variable $\eta$, the sea-surface height, |
---|
| 235 | is introduced which describes the shape of the air-sea interface. |
---|
[13463] | 236 | This variable is solution of a prognostic equation which is established by |
---|
| 237 | forming the vertical average of the kinematic surface condition (\autoref{eq:MB_w_bbc}): |
---|
[10414] | 238 | \begin{equation} |
---|
[13463] | 239 | \label{eq:MB_ssh} |
---|
[10501] | 240 | \pd[\eta]{t} = - D + P - E \quad \text{where} \quad D = \nabla \cdot \lt[ (H + \eta) \; \overline{U}_h \, \rt] |
---|
[707] | 241 | \end{equation} |
---|
[13463] | 242 | and using (\autoref{eq:MB_PE_hydrostatic}) the surface pressure is given by: |
---|
| 243 | $p_s = \rho \, g \, \eta$. |
---|
[707] | 244 | |
---|
[13463] | 245 | Allowing the air-sea interface to move introduces |
---|
| 246 | the \textbf{E}xternal \textbf{G}ravity \textbf{W}aves (EGWs) as |
---|
[10354] | 247 | a class of solution of the primitive equations. |
---|
[13463] | 248 | These waves are barotropic (\ie\ nearly independent of depth) and their phase speed is quite high. |
---|
[10354] | 249 | Their time scale is short with respect to the other processes described by the primitive equations. |
---|
[707] | 250 | |
---|
[10354] | 251 | Two choices can be made regarding the implementation of the free surface in the model, |
---|
[10501] | 252 | depending on the physical processes of interest. |
---|
[13463] | 253 | \begin{itemize} |
---|
| 254 | \item If one is interested in EGWs, in particular the tides and their interaction with |
---|
| 255 | the baroclinic structure of the ocean (internal waves) possibly in shallow seas, |
---|
| 256 | then a non linear free surface is the most appropriate. |
---|
| 257 | This means that no approximation is made in \autoref{eq:MB_ssh} and that |
---|
| 258 | the variation of the ocean volume is fully taken into account. |
---|
| 259 | Note that in order to study the fast time scales associated with EGWs it is necessary to |
---|
| 260 | minimize time filtering effects |
---|
| 261 | (use an explicit time scheme with very small time step, |
---|
| 262 | or a split-explicit scheme with reasonably small time step, |
---|
| 263 | see \autoref{subsec:DYN_spg_exp} or \autoref{subsec:DYN_spg_ts}). |
---|
| 264 | \item If one is not interested in EGWs but rather sees them as high frequency noise, |
---|
| 265 | it is possible to apply an explicit filter to slow down the fastest waves while |
---|
| 266 | not altering the slow barotropic Rossby waves. |
---|
| 267 | If further, an approximative conservation of heat and salt contents is sufficient for |
---|
| 268 | the problem solved, then it is sufficient to solve a linearized version of \autoref{eq:MB_ssh}, |
---|
| 269 | which still allows to take into account freshwater fluxes applied at the ocean surface |
---|
| 270 | \citep{roullet.madec_JGR00}. |
---|
| 271 | Nevertheless, with the linearization, an exact conservation of heat and salt contents is lost. |
---|
| 272 | \end{itemize} |
---|
| 273 | The filtering of EGWs in models with a free surface is usually |
---|
| 274 | a matter of discretisation of the temporal derivatives, |
---|
[11123] | 275 | using a split-explicit method \citep{killworth.webb.ea_JPO91, zhang.endoh_JGR92} or |
---|
[13463] | 276 | the implicit scheme \citep{dukowicz.smith_JGR94} or the addition of a filtering force in |
---|
| 277 | the momentum equation \citep{roullet.madec_JGR00}. |
---|
| 278 | With the present release, \NEMO\ offers the choice between an explicit free surface |
---|
| 279 | (see \autoref{subsec:DYN_spg_exp}) or a split-explicit scheme strongly inspired the one proposed by |
---|
| 280 | \citet{shchepetkin.mcwilliams_OM05} (see \autoref{subsec:DYN_spg_ts}). |
---|
[707] | 281 | |
---|
[13463] | 282 | %% ================================================================================================= |
---|
| 283 | \section{Curvilinear \textit{z}-coordinate system} |
---|
| 284 | \label{sec:MB_zco} |
---|
[707] | 285 | |
---|
[13463] | 286 | %% ================================================================================================= |
---|
[9393] | 287 | \subsection{Tensorial formalism} |
---|
[13463] | 288 | \label{subsec:MB_tensorial} |
---|
[707] | 289 | |
---|
[10354] | 290 | In many ocean circulation problems, the flow field has regions of enhanced dynamics |
---|
[13463] | 291 | (\ie\ surface layers, western boundary currents, equatorial currents, or ocean fronts). |
---|
[10354] | 292 | The representation of such dynamical processes can be improved by |
---|
| 293 | specifically increasing the model resolution in these regions. |
---|
| 294 | As well, it may be convenient to use a lateral boundary-following coordinate system to |
---|
| 295 | better represent coastal dynamics. |
---|
| 296 | Moreover, the common geographical coordinate system has a singular point at the North Pole that |
---|
| 297 | cannot be easily treated in a global model without filtering. |
---|
| 298 | A solution consists of introducing an appropriate coordinate transformation that |
---|
[11123] | 299 | shifts the singular point onto land \citep{madec.imbard_CD96, murray_JCP96}. |
---|
[13463] | 300 | As a consequence, |
---|
| 301 | it is important to solve the primitive equations in various curvilinear coordinate systems. |
---|
| 302 | An efficient way of introducing an appropriate coordinate transform can be found when |
---|
| 303 | using a tensorial formalism. |
---|
[10354] | 304 | This formalism is suited to any multidimensional curvilinear coordinate system. |
---|
[13463] | 305 | Ocean modellers mainly use three-dimensional orthogonal grids on the sphere |
---|
| 306 | (spherical earth approximation), with preservation of the local vertical. |
---|
| 307 | Here we give the simplified equations for this particular case. |
---|
| 308 | The general case is detailed by \citet{eiseman.stone_SR80} in |
---|
| 309 | their survey of the conservation laws of fluid dynamics. |
---|
[707] | 310 | |
---|
[10501] | 311 | Let $(i,j,k)$ be a set of orthogonal curvilinear coordinates on |
---|
[10354] | 312 | the sphere associated with the positively oriented orthogonal set of unit vectors |
---|
[10501] | 313 | $(i,j,k)$ linked to the earth such that |
---|
| 314 | $k$ is the local upward vector and $(i,j)$ are two vectors orthogonal to $k$, |
---|
[13463] | 315 | \ie\ along geopotential surfaces (\autoref{fig:MB_referential}). |
---|
[10354] | 316 | Let $(\lambda,\varphi,z)$ be the geographical coordinate system in which a position is defined by |
---|
| 317 | the latitude $\varphi(i,j)$, the longitude $\lambda(i,j)$ and |
---|
[10501] | 318 | the distance from the centre of the earth $a + z(k)$ where $a$ is the earth's radius and |
---|
[13463] | 319 | $z$ the altitude above a reference sea level (\autoref{fig:MB_referential}). |
---|
[10354] | 320 | The local deformation of the curvilinear coordinate system is given by $e_1$, $e_2$ and $e_3$, |
---|
| 321 | the three scale factors: |
---|
[10414] | 322 | \begin{equation} |
---|
[13463] | 323 | \label{eq:MB_scale_factors} |
---|
| 324 | e_1 = (a + z) \lt[ \lt( \pd[\lambda]{i} \cos \varphi \rt)^2 + \lt( \pd[\varphi]{i} \rt)^2 \rt]^{1/2} \quad e_2 = (a + z) \lt[ \lt( \pd[\lambda]{j} \cos \varphi \rt)^2 + \lt( \pd[\varphi]{j} \rt)^2 \rt]^{1/2} \quad e_3 = \lt( \pd[z]{k} \rt) |
---|
[10354] | 325 | \end{equation} |
---|
[707] | 326 | |
---|
[13463] | 327 | \begin{figure} |
---|
| 328 | \centering |
---|
| 329 | \includegraphics[width=0.33\textwidth]{MB_earth_referential} |
---|
| 330 | \caption[Geographical and curvilinear coordinate systems]{ |
---|
| 331 | the geographical coordinate system $(\lambda,\varphi,z)$ and the curvilinear |
---|
| 332 | coordinate system $(i,j,k)$.} |
---|
| 333 | \label{fig:MB_referential} |
---|
[10354] | 334 | \end{figure} |
---|
[707] | 335 | |
---|
[10501] | 336 | Since the ocean depth is far smaller than the earth's radius, $a + z$, can be replaced by $a$ in |
---|
[13463] | 337 | (\autoref{eq:MB_scale_factors}) (thin-shell approximation). |
---|
[10354] | 338 | The resulting horizontal scale factors $e_1$, $e_2$ are independent of $k$ while |
---|
[10501] | 339 | the vertical scale factor is a single function of $k$ as $k$ is parallel to $z$. |
---|
[10354] | 340 | The scalar and vector operators that appear in the primitive equations |
---|
[13463] | 341 | (\autoref{eq:MB_PE_dyn} to \autoref{eq:MB_PE_eos}) can then be written in the tensorial form, |
---|
[10354] | 342 | invariant in any orthogonal horizontal curvilinear coordinate system transformation: |
---|
[10414] | 343 | \begin{subequations} |
---|
[13463] | 344 | % \label{eq:MB_discrete_operators} |
---|
| 345 | \begin{align} |
---|
| 346 | \label{eq:MB_grad} |
---|
| 347 | \nabla q &= \frac{1}{e_1} \pd[q]{i} \; \vect i + \frac{1}{e_2} \pd[q]{j} \; \vect j + \frac{1}{e_3} \pd[q]{k} \; \vect k \\ |
---|
| 348 | \label{eq:MB_div} |
---|
| 349 | \nabla \cdot \vect A &= \frac{1}{e_1 \; e_2} \lt[ \pd[(e_2 \; a_1)]{\partial i} + \pd[(e_1 \; a_2)]{j} \rt] + \frac{1}{e_3} \lt[ \pd[a_3]{k} \rt] \\ |
---|
| 350 | \label{eq:MB_curl} |
---|
| 351 | \nabla \times \vect{A} &= \lt[ \frac{1}{e_2} \pd[a_3]{j} - \frac{1}{e_3} \pd[a_2]{k} \rt] \vect i + \lt[ \frac{1}{e_3} \pd[a_1]{k} - \frac{1}{e_1} \pd[a_3]{i} \rt] \vect j + \frac{1}{e_1 e_2} \lt[ \pd[(e_2 a_2)]{i} - \pd[(e_1 a_1)]{j} \rt] \vect k \\ |
---|
| 352 | \label{eq:MB_lap} |
---|
| 353 | \Delta q &= \nabla \cdot (\nabla q) \\ |
---|
| 354 | \label{eq:MB_lap_vector} |
---|
| 355 | \Delta \vect A &= \nabla (\nabla \cdot \vect A) - \nabla \times (\nabla \times \vect A) |
---|
| 356 | \end{align} |
---|
[707] | 357 | \end{subequations} |
---|
[13463] | 358 | where $q$ is a scalar quantity and |
---|
| 359 | $\vect A = (a_1,a_2,a_3)$ a vector in the $(i,j,k)$ coordinates system. |
---|
[707] | 360 | |
---|
[13463] | 361 | %% ================================================================================================= |
---|
[9393] | 362 | \subsection{Continuous model equations} |
---|
[13463] | 363 | \label{subsec:MB_zco_Eq} |
---|
[707] | 364 | |
---|
[10354] | 365 | In order to express the Primitive Equations in tensorial formalism, |
---|
[13463] | 366 | it is necessary to compute the horizontal component of the non-linear and viscous terms of |
---|
| 367 | the equation using \autoref{eq:MB_grad}) to \autoref{eq:MB_lap_vector}. |
---|
| 368 | Let us set $\vect U = (u,v,w) = \vect U_h + w \; \vect k $, |
---|
| 369 | the velocity in the $(i,j,k)$ coordinates system, |
---|
| 370 | and define the relative vorticity $\zeta$ and the divergence of the horizontal velocity field $\chi$, |
---|
| 371 | by: |
---|
[10501] | 372 | \begin{gather} |
---|
[13463] | 373 | \label{eq:MB_curl_Uh} |
---|
[10501] | 374 | \zeta = \frac{1}{e_1 e_2} \lt[ \pd[(e_2 \, v)]{i} - \pd[(e_1 \, u)]{j} \rt] \\ |
---|
[13463] | 375 | \label{eq:MB_div_Uh} |
---|
[10501] | 376 | \chi = \frac{1}{e_1 e_2} \lt[ \pd[(e_2 \, u)]{i} + \pd[(e_1 \, v)]{j} \rt] |
---|
| 377 | \end{gather} |
---|
[707] | 378 | |
---|
[11335] | 379 | Using again the fact that the horizontal scale factors $e_1$ and $e_2$ are independent of $k$ and that |
---|
[10354] | 380 | $e_3$ is a function of the single variable $k$, |
---|
[13463] | 381 | $NLT$ the nonlinear term of \autoref{eq:MB_PE_dyn} can be transformed as follows: |
---|
| 382 | \begin{align*} |
---|
| 383 | NLT &= \lt[ (\nabla \times {\vect U}) \times {\vect U} + \frac{1}{2} \nabla \lt( {\vect U}^2 \rt) \rt]_h \\ |
---|
| 384 | &= \lt( |
---|
| 385 | \begin{array}{*{20}c} |
---|
| 386 | \lt[ \frac{1}{e_3} \pd[u]{k} - \frac{1}{e_1} \pd[w]{i} \rt] w - \zeta \; v \\ |
---|
| 387 | \zeta \; u - \lt[ \frac{1}{e_2} \pd[w]{j} - \frac{1}{e_3} \pd[v]{k} \rt] \ w |
---|
| 388 | \end{array} |
---|
| 389 | \rt) |
---|
| 390 | + \frac{1}{2} \lt( |
---|
| 391 | \begin{array}{*{20}c} |
---|
| 392 | \frac{1}{e_1} \pd[(u^2 + v^2 + w^2)]{i} \\ |
---|
| 393 | \frac{1}{e_2} \pd[(u^2 + v^2 + w^2)]{j} |
---|
| 394 | \end{array} |
---|
| 395 | \rt) \\ |
---|
| 396 | &= \lt( |
---|
| 397 | \begin{array}{*{20}c} |
---|
| 398 | -\zeta \; v \\ |
---|
| 399 | \zeta \; u |
---|
| 400 | \end{array} |
---|
| 401 | \rt) |
---|
| 402 | + \frac{1}{2} \lt( |
---|
| 403 | \begin{array}{*{20}c} |
---|
| 404 | \frac{1}{e_1} \pd[(u^2 + v^2)]{i} \\ |
---|
| 405 | \frac{1}{e_2} \pd[(u^2 + v^2)]{j} |
---|
| 406 | \end{array} |
---|
| 407 | \rt) |
---|
| 408 | + \frac{1}{e_3} \lt( |
---|
| 409 | \begin{array}{*{20}c} |
---|
| 410 | w \; \pd[u]{k} \\ |
---|
| 411 | w \; \pd[v]{k} |
---|
| 412 | \end{array} |
---|
| 413 | \rt) |
---|
| 414 | - \lt( |
---|
| 415 | \begin{array}{*{20}c} |
---|
| 416 | \frac{w}{e_1} \pd[w]{i} - \frac{1}{2 e_1} \pd[w^2]{i} \\ |
---|
| 417 | \frac{w}{e_2} \pd[w]{j} - \frac{1}{2 e_2} \pd[w^2]{j} |
---|
| 418 | \end{array} |
---|
| 419 | \rt) |
---|
| 420 | \end{align*} |
---|
| 421 | The last term of the right hand side is obviously zero, |
---|
| 422 | and thus the \textbf{N}on\textbf{L}inear \textbf{T}erm ($NLT$) of \autoref{eq:MB_PE_dyn} is written in |
---|
| 423 | the $(i,j,k)$ coordinate system: |
---|
[10414] | 424 | \begin{equation} |
---|
[13463] | 425 | \label{eq:MB_vector_form} |
---|
| 426 | NLT = \zeta \; \vect k \times \vect U_h + \frac{1}{2} \nabla_h \lt( \vect U_h^2 \rt) + \frac{1}{e_3} w \pd[\vect U_h]{k} |
---|
[707] | 427 | \end{equation} |
---|
| 428 | |
---|
[10354] | 429 | This is the so-called \textit{vector invariant form} of the momentum advection term. |
---|
| 430 | For some purposes, it can be advantageous to write this term in the so-called flux form, |
---|
[13463] | 431 | \ie\ to write it as the divergence of fluxes. |
---|
| 432 | For example, |
---|
| 433 | the first component of \autoref{eq:MB_vector_form} (the $i$-component) is transformed as follows: |
---|
| 434 | \begin{alignat*}{3} |
---|
[10501] | 435 | &NLT_i &= &- \zeta \; v + \frac{1}{2 \; e_1} \pd[ (u^2 + v^2) ]{i} + \frac{1}{e_3} w \ \pd[u]{k} \\ |
---|
[13463] | 436 | & &= &\frac{1}{e_1 \; e_2} \lt( -v \pd[(e_2 \, v)]{i} + v \pd[(e_1 \, u)]{j} \rt) + \frac{1}{e_1 e_2} \lt( e_2 \; u \pd[u]{i} + e_2 \; v \pd[v]{i} \rt) + \frac{1}{e_3} \lt( w \; \pd[u]{k} \rt) \\ |
---|
| 437 | & &= &\frac{1}{e_1 \; e_2} \lt[ - \lt( v^2 \pd[e_2]{i} + e_2 \, v \pd[v]{i} \rt) + \lt( \pd[ \lt( e_1 \, u \, v \rt)]{j} - e_1 \, u \pd[v]{j} \rt) \rt. \lt. + \lt( \pd[ \lt( e_2 \, u \, u \rt)]{i} - u \pd[ \lt( e_2 u \rt)]{i} \rt) + e_2 v \pd[v]{i} \rt] \\ |
---|
[10501] | 438 | & & &+ \frac{1}{e_3} \lt( \pd[(w \, u)]{k} - u \pd[w]{k} \rt) \\ |
---|
[13463] | 439 | & &= &\frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, u)]{i} + \pd[(e_1 \, u \, v)]{j} \rt) + \frac{1}{e_3} \pd[(w \, u)]{k} + \frac{1}{e_1 e_2} \lt[ - u \lt( \pd[(e_1 v)]{j} - v \, \pd[e_1]{j} \rt) - u \pd[(e_2 u)]{i} \rt] - \frac{1}{e_3} \pd[w]{k} u \\ |
---|
[10501] | 440 | & & &+ \frac{1}{e_1 e_2} \lt( - v^2 \pd[e_2]{i} \rt) \\ |
---|
[13463] | 441 | & &= &\nabla \cdot (\vect U \, u) - (\nabla \cdot \vect U) \ u + \frac{1}{e_1 e_2} \lt( -v^2 \pd[e_2]{i} + u v \, \pd[e_1]{j} \rt) \\ |
---|
| 442 | \shortintertext{as $\nabla \cdot {\vect U} \; = 0$ (incompressibility) it becomes:} |
---|
[10501] | 443 | & &= &\, \nabla \cdot (\vect U \, u) + \frac{1}{e_1 e_2} \lt( v \; \pd[e_2]{i} - u \; \pd[e_1]{j} \rt) (-v) |
---|
| 444 | \end{alignat*} |
---|
[707] | 445 | |
---|
[817] | 446 | The flux form of the momentum advection term is therefore given by: |
---|
[10501] | 447 | \begin{equation} |
---|
[13463] | 448 | \label{eq:MB_flux_form} |
---|
| 449 | NLT = \nabla \cdot \lt( |
---|
| 450 | \begin{array}{*{20}c} |
---|
| 451 | \vect U \, u \\ |
---|
| 452 | \vect U \, v |
---|
| 453 | \end{array} |
---|
| 454 | \rt) |
---|
| 455 | + \frac{1}{e_1 e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \vect k \times \vect U_h |
---|
[10501] | 456 | \end{equation} |
---|
[707] | 457 | |
---|
[10354] | 458 | The flux form has two terms, |
---|
[13463] | 459 | the first one is expressed as the divergence of momentum fluxes |
---|
| 460 | (hence the flux form name given to this formulation) and |
---|
| 461 | the second one is due to the curvilinear nature of the coordinate system used. |
---|
| 462 | The latter is called the \textit{metric} term and can be viewed as |
---|
| 463 | a modification of the Coriolis parameter: |
---|
[10414] | 464 | \[ |
---|
[13463] | 465 | % \label{eq:MB_cor+metric} |
---|
[10501] | 466 | f \to f + \frac{1}{e_1 e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) |
---|
[10414] | 467 | \] |
---|
[707] | 468 | |
---|
[10354] | 469 | Note that in the case of geographical coordinate, |
---|
[13463] | 470 | \ie\ when $(i,j) \to (\lambda,\varphi)$ and $(e_1,e_2) \to (a \, \cos \varphi,a)$, |
---|
[10501] | 471 | we recover the commonly used modification of the Coriolis parameter $f \to f + (u / a) \tan \varphi$. |
---|
[707] | 472 | |
---|
[10354] | 473 | To sum up, the curvilinear $z$-coordinate equations solved by the ocean model can be written in |
---|
| 474 | the following tensorial formalism: |
---|
[2282] | 475 | |
---|
[13463] | 476 | \begin{description} |
---|
| 477 | \item [Vector invariant form of the momentum equations] |
---|
[10501] | 478 | \begin{equation} |
---|
[13463] | 479 | \label{eq:MB_dyn_vect} |
---|
| 480 | \begin{gathered} |
---|
| 481 | % \label{eq:MB_dyn_vect_u} |
---|
| 482 | \pd[u]{t} = + (\zeta + f) \, v - \frac{1}{2 e_1} \pd[]{i} (u^2 + v^2) - \frac{1}{e_3} w \pd[u]{k} - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) + D_u^{\vect U} + F_u^{\vect U} \\ |
---|
| 483 | \pd[v]{t} = - (\zeta + f) \, u - \frac{1}{2 e_2} \pd[]{j} (u^2 + v^2) - \frac{1}{e_3} w \pd[v]{k} - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) + D_v^{\vect U} + F_v^{\vect U} |
---|
| 484 | \end{gathered} |
---|
[10501] | 485 | \end{equation} |
---|
[13463] | 486 | \item [Flux form of the momentum equations] |
---|
| 487 | % \label{eq:MB_dyn_flux} |
---|
| 488 | \begin{alignat*}{2} |
---|
| 489 | % \label{eq:MB_dyn_flux_u} |
---|
| 490 | \pd[u]{t} = &+ \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, v - \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, u)]{i} + \pd[(e_1 \, v \, u)]{j} \rt) - \frac{1}{e_3} \pd[(w \, u)]{k} \\ |
---|
| 491 | &- \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) + D_u^{\vect U} + F_u^{\vect U} |
---|
| 492 | \end{alignat*} |
---|
| 493 | \begin{alignat*}{2} |
---|
| 494 | % \label{eq:MB_dyn_flux_v} |
---|
| 495 | \pd[v]{t} = &- \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, u - \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, v)]{i} + \pd[(e_1 \, v \, v)]{j} \rt) - \frac{1}{e_3} \pd[(w \, v)]{k} \\ |
---|
| 496 | &- \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) + D_v^{\vect U} + F_v^{\vect U} |
---|
| 497 | \end{alignat*} |
---|
| 498 | where $\zeta$, the relative vorticity, is given by \autoref{eq:MB_curl_Uh} and |
---|
| 499 | $p_s$, the surface pressure, is given by: |
---|
[10501] | 500 | \[ |
---|
[13463] | 501 | % \label{eq:MB_spg} |
---|
[10501] | 502 | p_s = \rho \,g \, \eta |
---|
| 503 | \] |
---|
[13463] | 504 | and $\eta$ is the solution of \autoref{eq:MB_ssh}. |
---|
[707] | 505 | |
---|
[10501] | 506 | The vertical velocity and the hydrostatic pressure are diagnosed from the following equations: |
---|
| 507 | \[ |
---|
[13463] | 508 | % \label{eq:MB_w_diag} |
---|
| 509 | \pd[w]{k} = - \chi \; e_3 \qquad |
---|
| 510 | % \label{eq:MB_hp_diag} |
---|
[10501] | 511 | \pd[p_h]{k} = - \rho \; g \; e_3 |
---|
| 512 | \] |
---|
[13463] | 513 | where the divergence of the horizontal velocity, $\chi$ is given by \autoref{eq:MB_div_Uh}. |
---|
| 514 | \item [Tracer equations] |
---|
| 515 | \begin{gather*} |
---|
| 516 | \pd[T]{t} = - \frac{1}{e_1 e_2} \lt[ \pd[(e_2 T \, u)]{i} + \pd[(e_1 T \, v)]{j} \rt] - \frac{1}{e_3} \pd[(T \, w)]{k} + D^T + F^T \\ |
---|
| 517 | \pd[S]{t} = - \frac{1}{e_1 e_2} \lt[ \pd[(e_2 S \, u)]{i} + \pd[(e_1 S \, v)]{j} \rt] - \frac{1}{e_3} \pd[(S \, w)]{k} + D^S + F^S \\ |
---|
| 518 | \rho = \rho \big( T,S,z(k) \big) |
---|
| 519 | \end{gather*} |
---|
| 520 | \end{description} |
---|
[11335] | 521 | |
---|
[13463] | 522 | The expression of $\vect D^{U}$, $D^{S}$ and $D^{T}$ depends on |
---|
| 523 | the subgrid scale parameterisation used. |
---|
| 524 | It will be defined in \autoref{eq:MB_zdf}. |
---|
[10501] | 525 | The nature and formulation of $\vect F^{\vect U}$, $F^T$ and $F^S$, the surface forcing terms, |
---|
[10354] | 526 | are discussed in \autoref{chap:SBC}. |
---|
[707] | 527 | |
---|
[13463] | 528 | %% ================================================================================================= |
---|
[9393] | 529 | \section{Curvilinear generalised vertical coordinate system} |
---|
[13463] | 530 | \label{sec:MB_gco} |
---|
[817] | 531 | |
---|
[10354] | 532 | The ocean domain presents a huge diversity of situation in the vertical. |
---|
| 533 | First the ocean surface is a time dependent surface (moving surface). |
---|
| 534 | Second the ocean floor depends on the geographical position, |
---|
| 535 | varying from more than 6,000 meters in abyssal trenches to zero at the coast. |
---|
[10501] | 536 | Last but not least, the ocean stratification exerts a strong barrier to vertical motions and mixing. |
---|
[13463] | 537 | Therefore, in order to represent the ocean with respect to the first point |
---|
| 538 | a space and time dependent vertical coordinate that follows the variation of the sea surface height |
---|
| 539 | \eg\ an \zstar-coordinate; |
---|
| 540 | for the second point, |
---|
| 541 | a space variation to fit the change of bottom topography |
---|
| 542 | \eg\ a terrain-following or $\sigma$-coordinate; |
---|
| 543 | and for the third point, |
---|
| 544 | one will be tempted to use a space and time dependent coordinate that follows the isopycnal surfaces, |
---|
| 545 | \eg\ an isopycnic coordinate. |
---|
[2282] | 546 | |
---|
[13463] | 547 | In order to satisfy two or more constraints one can even be tempted to mixed these coordinate systems, |
---|
| 548 | as in HYCOM (mixture of $z$-coordinate at the surface, isopycnic coordinate in the ocean interior and |
---|
| 549 | $\sigma$ at the ocean bottom) \citep{chassignet.smith.ea_JPO03} or |
---|
| 550 | OPA (mixture of $z$-coordinate in vicinity the surface and steep topography areas and |
---|
| 551 | $\sigma$-coordinate elsewhere) \citep{madec.delecluse.ea_JPO96} among others. |
---|
[2282] | 552 | |
---|
[10354] | 553 | In fact one is totally free to choose any space and time vertical coordinate by |
---|
| 554 | introducing an arbitrary vertical coordinate : |
---|
[10414] | 555 | \begin{equation} |
---|
[13463] | 556 | \label{eq:MB_s} |
---|
[10501] | 557 | s = s(i,j,k,t) |
---|
[2282] | 558 | \end{equation} |
---|
[13463] | 559 | with the restriction that the above equation gives a single-valued monotonic relationship between |
---|
| 560 | $s$ and $k$, when $i$, $j$ and $t$ are held fixed. |
---|
| 561 | \autoref{eq:MB_s} is a transformation from |
---|
| 562 | the $(i,j,k,t)$ coordinate system with independent variables into |
---|
[10354] | 563 | the $(i,j,s,t)$ generalised coordinate system with $s$ depending on the other three variables through |
---|
[13463] | 564 | \autoref{eq:MB_s}. |
---|
[11123] | 565 | This so-called \textit{generalised vertical coordinate} \citep{kasahara_MWR74} is in fact |
---|
[13463] | 566 | an \textbf{A}rbitrary \textbf{L}agrangian--\textbf{E}ulerian (ALE) coordinate. |
---|
| 567 | Indeed, one has a great deal of freedom in the choice of expression for $s$. |
---|
| 568 | The choice determines which part of the vertical velocity (defined from a fixed referential) |
---|
| 569 | will cross the levels (Eulerian part) and which part will be used to move them (Lagrangian part). |
---|
| 570 | The coordinate is also sometimes referenced as an adaptive coordinate |
---|
| 571 | \citep{hofmeister.burchard.ea_OM10}, since the coordinate system is adapted in |
---|
| 572 | the course of the simulation. |
---|
[10354] | 573 | Its most often used implementation is via an ALE algorithm, |
---|
| 574 | in which a pure lagrangian step is followed by regridding and remapping steps, |
---|
[11335] | 575 | the latter step implicitly embedding the vertical advection |
---|
[11123] | 576 | \citep{hirt.amsden.ea_JCP74, chassignet.smith.ea_JPO03, white.adcroft.ea_JCP09}. |
---|
| 577 | Here we follow the \citep{kasahara_MWR74} strategy: |
---|
[11335] | 578 | a regridding step (an update of the vertical coordinate) followed by an Eulerian step with |
---|
[10354] | 579 | an explicit computation of vertical advection relative to the moving s-surfaces. |
---|
[2282] | 580 | |
---|
[13463] | 581 | \cmtgm{A key point here is that the $s$-coordinate depends on $(i,j)$ |
---|
| 582 | ==> horizontal pressure gradient...} |
---|
[11335] | 583 | The generalized vertical coordinates used in ocean modelling are not orthogonal, |
---|
[10354] | 584 | which contrasts with many other applications in mathematical physics. |
---|
| 585 | Hence, it is useful to keep in mind the following properties that may seem odd on initial encounter. |
---|
[2282] | 586 | |
---|
[10354] | 587 | The horizontal velocity in ocean models measures motions in the horizontal plane, |
---|
| 588 | perpendicular to the local gravitational field. |
---|
[13463] | 589 | That is, horizontal velocity is mathematically the same regardless of the vertical coordinate, |
---|
| 590 | be it geopotential, isopycnal, pressure, or terrain following. |
---|
[10354] | 591 | The key motivation for maintaining the same horizontal velocity component is that |
---|
| 592 | the hydrostatic and geostrophic balances are dominant in the large-scale ocean. |
---|
[13463] | 593 | Use of an alternative quasi -horizontal velocity, |
---|
| 594 | for example one oriented parallel to the generalized surface, |
---|
[10354] | 595 | would lead to unacceptable numerical errors. |
---|
[10501] | 596 | Correspondingly, the vertical direction is anti -parallel to the gravitational force in |
---|
[10354] | 597 | all of the coordinate systems. |
---|
[10501] | 598 | We do not choose the alternative of a quasi -vertical direction oriented normal to |
---|
| 599 | the surface of a constant generalized vertical coordinate. |
---|
[2282] | 600 | |
---|
[13463] | 601 | It is the method used to measure transport across the generalized vertical coordinate surfaces which |
---|
| 602 | differs between the vertical coordinate choices. |
---|
| 603 | That is, |
---|
| 604 | computation of the dia-surface velocity component represents the fundamental distinction between |
---|
[10354] | 605 | the various coordinates. |
---|
[13463] | 606 | In some models, such as geopotential, pressure, and terrain following, |
---|
| 607 | this transport is typically diagnosed from volume or mass conservation. |
---|
| 608 | In other models, such as isopycnal layered models, |
---|
| 609 | this transport is prescribed based on assumptions about the physical processes producing |
---|
| 610 | a flux across the layer interfaces. |
---|
[2282] | 611 | |
---|
[10354] | 612 | In this section we first establish the PE in the generalised vertical $s$-coordinate, |
---|
[10501] | 613 | then we discuss the particular cases available in \NEMO, namely $z$, \zstar, $s$, and \ztilde. |
---|
[2282] | 614 | %} |
---|
| 615 | |
---|
[13463] | 616 | %% ================================================================================================= |
---|
[10501] | 617 | \subsection{\textit{S}-coordinate formulation} |
---|
[2282] | 618 | |
---|
[13463] | 619 | Starting from the set of equations established in \autoref{sec:MB_zco} for |
---|
| 620 | the special case $k = z$ and thus $e_3 = 1$, |
---|
| 621 | we introduce an arbitrary vertical coordinate $s = s(i,j,k,t)$, |
---|
[10501] | 622 | which includes $z$-, \zstar- and $\sigma$-coordinates as special cases |
---|
| 623 | ($s = z$, $s = \zstar$, and $s = \sigma = z / H$ or $ = z / \lt( H + \eta \rt)$, resp.). |
---|
[13463] | 624 | A formal derivation of the transformed equations is given in \autoref{apdx:SCOORD}. |
---|
| 625 | Let us define the vertical scale factor by $e_3 = \partial_s z$ |
---|
| 626 | ($e_3$ is now a function of $(i,j,k,t)$ ), |
---|
[10501] | 627 | and the slopes in the $(i,j)$ directions between $s$- and $z$-surfaces by: |
---|
[10414] | 628 | \begin{equation} |
---|
[13463] | 629 | \label{eq:MB_sco_slope} |
---|
[10501] | 630 | \sigma_1 = \frac{1}{e_1} \; \lt. \pd[z]{i} \rt|_s \quad \text{and} \quad |
---|
| 631 | \sigma_2 = \frac{1}{e_2} \; \lt. \pd[z]{j} \rt|_s |
---|
[2282] | 632 | \end{equation} |
---|
[13463] | 633 | We also introduce $\omega$, a dia-surface velocity component, |
---|
| 634 | defined as the velocity relative to the moving $s$-surfaces and normal to them: |
---|
[10414] | 635 | \[ |
---|
[13463] | 636 | % \label{eq:MB_sco_w} |
---|
[11335] | 637 | \omega = w - \, \lt. \pd[z]{t} \rt|_s - \sigma_1 \, u - \sigma_2 \, v |
---|
[10414] | 638 | \] |
---|
[2282] | 639 | |
---|
[13463] | 640 | The equations solved by the ocean model \autoref{eq:MB_PE} in $s$-coordinate can be written as follows |
---|
| 641 | (see \autoref{sec:SCOORD_momentum}): |
---|
[2282] | 642 | |
---|
[13463] | 643 | \begin{description} |
---|
| 644 | \item [Vector invariant form of the momentum equation] |
---|
| 645 | \begin{gather*} |
---|
| 646 | % \label{eq:MB_sco_u_vector} |
---|
| 647 | \pd[u]{t} = + (\zeta + f) \, v - \frac{1}{2 \, e_1} \pd[]{i} (u^2 + v^2) - \frac{1}{e_3} \omega \pd[u]{k} - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) - g \frac{\rho}{\rho_o} \sigma_1 + D_u^{\vect U} + F_u^{\vect U} \\ |
---|
| 648 | % \label{eq:MB_sco_v_vector} |
---|
| 649 | \pd[v]{t} = - (\zeta + f) \, u - \frac{1}{2 \, e_2} \pd[]{j} (u^2 + v^2) - \frac{1}{e_3} \omega \pd[v]{k} - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) - g \frac{\rho}{\rho_o} \sigma_2 + D_v^{\vect U} + F_v^{\vect U} |
---|
| 650 | \end{gather*} |
---|
| 651 | \item [Flux form of the momentum equation] |
---|
| 652 | \begin{alignat*}{2} |
---|
| 653 | % \label{eq:MB_sco_u_flux} |
---|
| 654 | \frac{1}{e_3} \pd[(e_3 \, u)]{t} = &+ \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, v - \frac{1}{e_1 \; e_2 \; e_3} \lt( \pd[(e_2 \, e_3 \, u \, u)]{i} + \pd[(e_1 \, e_3 \, v \, u)]{j} \rt) - \frac{1}{e_3} \pd[(\omega \, u)]{k} \\ |
---|
| 655 | &- \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) - g \frac{\rho}{\rho_o} \sigma_1 + D_u^{\vect U} + F_u^{\vect U} |
---|
| 656 | \end{alignat*} |
---|
| 657 | \begin{alignat*}{2} |
---|
| 658 | % \label{eq:MB_sco_v_flux} |
---|
| 659 | \frac{1}{e_3} \pd[(e_3 \, v)]{t} = &- \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, u - \frac{1}{e_1 \; e_2 \; e_3} \lt( \pd[( e_2 \; e_3 \, u \, v)]{i} + \pd[(e_1 \; e_3 \, v \, v)]{j} \rt) - \frac{1}{e_3} \pd[(\omega \, v)]{k} \\ |
---|
| 660 | &- \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) - g \frac{\rho}{\rho_o}\sigma_2 + D_v^{\vect U} + F_v^{\vect U} |
---|
| 661 | \end{alignat*} |
---|
[10501] | 662 | where the relative vorticity, $\zeta$, the surface pressure gradient, |
---|
| 663 | and the hydrostatic pressure have the same expressions as in $z$-coordinates although |
---|
| 664 | they do not represent exactly the same quantities. |
---|
[13463] | 665 | $\omega$ is provided by the continuity equation (see \autoref{apdx:SCOORD}): |
---|
[10501] | 666 | \[ |
---|
[13463] | 667 | % \label{eq:MB_sco_continuity} |
---|
[10501] | 668 | \pd[e_3]{t} + e_3 \; \chi + \pd[\omega]{s} = 0 \quad \text{with} \quad |
---|
| 669 | \chi = \frac{1}{e_1 e_2 e_3} \lt( \pd[(e_2 e_3 \, u)]{i} + \pd[(e_1 e_3 \, v)]{j} \rt) |
---|
| 670 | \] |
---|
[13463] | 671 | \item [Tracer equations] |
---|
| 672 | \begin{gather*} |
---|
| 673 | % \label{eq:MB_sco_t} |
---|
| 674 | \frac{1}{e_3} \pd[(e_3 \, T)]{t} = - \frac{1}{e_1 e_2 e_3} \lt( \pd[(e_2 e_3 \, u \, T)]{i} + \pd[(e_1 e_3 \, v \, T)]{j} \rt) - \frac{1}{e_3} \pd[(T \, \omega)]{k} + D^T + F^S \\ |
---|
| 675 | % \label{eq:MB_sco_s} |
---|
| 676 | \frac{1}{e_3} \pd[(e_3 \, S)]{t} = - \frac{1}{e_1 e_2 e_3} \lt( \pd[(e_2 e_3 \, u \, S)]{i} + \pd[(e_1 e_3 \, v \, S)]{j} \rt) - \frac{1}{e_3} \pd[(S \, \omega)]{k} + D^S + F^S |
---|
| 677 | \end{gather*} |
---|
| 678 | \end{description} |
---|
[10354] | 679 | The equation of state has the same expression as in $z$-coordinate, |
---|
| 680 | and similar expressions are used for mixing and forcing terms. |
---|
[2282] | 681 | |
---|
[13463] | 682 | \cmtgm{ |
---|
[10501] | 683 | \colorbox{yellow}{ to be updated $= = >$} |
---|
| 684 | Add a few works on z and zps and s and underlies the differences between all of them |
---|
| 685 | \colorbox{yellow}{$< = =$ end update} |
---|
| 686 | } |
---|
[2282] | 687 | |
---|
[13463] | 688 | %% ================================================================================================= |
---|
[10501] | 689 | \subsection{Curvilinear \zstar-coordinate system} |
---|
[13463] | 690 | \label{subsec:MB_zco_star} |
---|
[2282] | 691 | |
---|
[13463] | 692 | \begin{figure} |
---|
| 693 | \centering |
---|
| 694 | \includegraphics[width=0.66\textwidth]{MB_z_zstar} |
---|
| 695 | \caption[Curvilinear z-coordinate systems (\{non-\}linear free-surface cases and re-scaled \zstar)]{ |
---|
| 696 | \begin{enumerate*}[label=(\textit{\alph*})] |
---|
| 697 | \item $z$-coordinate in linear free-surface case; |
---|
| 698 | \item $z$-coordinate in non-linear free surface case; |
---|
| 699 | \item re-scaled height coordinate |
---|
[11123] | 700 | (become popular as the \zstar-coordinate \citep{adcroft.campin_OM04}). |
---|
[13463] | 701 | \end{enumerate*} |
---|
| 702 | } |
---|
| 703 | \label{fig:MB_z_zstar} |
---|
[10354] | 704 | \end{figure} |
---|
[817] | 705 | |
---|
[13463] | 706 | In this case, the free surface equation is nonlinear, |
---|
| 707 | and the variations of volume are fully taken into account. |
---|
| 708 | These coordinates systems is presented in a report \citep{levier.treguier.ea_rpt07} available on |
---|
| 709 | the \NEMO\ web site. |
---|
[817] | 710 | |
---|
[13463] | 711 | The \zstar\ coordinate approach is an unapproximated, non-linear free surface implementation which |
---|
| 712 | allows one to deal with large amplitude free-surface variations relative to the vertical resolution |
---|
| 713 | \citep{adcroft.campin_OM04}. |
---|
| 714 | In the \zstar\ formulation, the variation of the column thickness due to sea-surface undulations is |
---|
| 715 | not concentrated in the surface level, as in the $z$-coordinate formulation, |
---|
| 716 | but is equally distributed over the full water column. |
---|
[10354] | 717 | Thus vertical levels naturally follow sea-surface variations, with a linear attenuation with depth, |
---|
[13463] | 718 | as illustrated by \autoref{fig:MB_z_zstar}. |
---|
| 719 | Note that with a flat bottom, such as in \autoref{fig:MB_z_zstar}, |
---|
| 720 | the bottom-following $z$ coordinate and \zstar\ are equivalent. |
---|
[10501] | 721 | The definition and modified oceanic equations for the rescaled vertical coordinate \zstar, |
---|
[10354] | 722 | including the treatment of fresh-water flux at the surface, are detailed in Adcroft and Campin (2004). |
---|
| 723 | The major points are summarized here. |
---|
[13463] | 724 | The position (\zstar) and vertical discretization ($\delta \zstar$) are expressed as: |
---|
[10414] | 725 | \[ |
---|
[13463] | 726 | % \label{eq:MB_z-star} |
---|
| 727 | H + \zstar = (H + z) / r \quad \text{and} \quad \delta \zstar = \delta z / r \quad \text{with} \quad r = \frac{H + \eta}{H} |
---|
[10501] | 728 | \] |
---|
[11335] | 729 | Simple re-organisation of the above expressions gives |
---|
| 730 | \[ |
---|
[13463] | 731 | % \label{eq:MB_zstar_2} |
---|
| 732 | \zstar = H \lt( \frac{z - \eta}{H + \eta} \rt) |
---|
[11335] | 733 | \] |
---|
[10442] | 734 | Since the vertical displacement of the free surface is incorporated in the vertical coordinate \zstar, |
---|
[13463] | 735 | the upper and lower boundaries are at fixed \zstar\ position, |
---|
[10501] | 736 | $\zstar = 0$ and $\zstar = -H$ respectively. |
---|
[10354] | 737 | Also the divergence of the flow field is no longer zero as shown by the continuity equation: |
---|
[10501] | 738 | \[ |
---|
[13463] | 739 | \pd[r]{t} = \nabla_{\zstar} \cdot \lt( r \; \vect U_h \rt) + \pd[r \; w^*]{\zstar} = 0 |
---|
[10501] | 740 | \] |
---|
[13463] | 741 | This \zstar\ coordinate is closely related to the $\eta$ coordinate used in many atmospheric models |
---|
| 742 | (see Black (1994) for a review of $\eta$ coordinate atmospheric models). |
---|
[10354] | 743 | It was originally used in ocean models by Stacey et al. (1995) for studies of tides next to shelves, |
---|
| 744 | and it has been recently promoted by Adcroft and Campin (2004) for global climate modelling. |
---|
[2282] | 745 | |
---|
[13463] | 746 | The surfaces of constant \zstar\ are quasi-horizontal. |
---|
| 747 | Indeed, the \zstar\ coordinate reduces to $z$ when $\eta$ is zero. |
---|
[10354] | 748 | In general, when noting the large differences between |
---|
| 749 | undulations of the bottom topography versus undulations in the surface height, |
---|
[13463] | 750 | it is clear that surfaces constant \zstar\ are very similar to the depth surfaces. |
---|
[10354] | 751 | These properties greatly reduce difficulties of computing the horizontal pressure gradient relative to |
---|
[13463] | 752 | terrain following sigma models discussed in \autoref{subsec:MB_sco}. |
---|
[11335] | 753 | Additionally, since $\zstar = z$ when $\eta = 0$, |
---|
[13463] | 754 | no flow is spontaneously generated in an unforced ocean starting from rest, |
---|
| 755 | regardless the bottom topography. |
---|
| 756 | This behaviour is in contrast to the case with ``s''-models, |
---|
| 757 | where pressure gradient errors in the presence of nontrivial topographic variations can generate |
---|
| 758 | nontrivial spontaneous flow from a resting state, |
---|
[10354] | 759 | depending on the sophistication of the pressure gradient solver. |
---|
[13463] | 760 | The quasi-horizontal nature of the coordinate surfaces also facilitates the implementation of |
---|
| 761 | neutral physics parameterizations in \zstar\ models using the same techniques as in $z$-models |
---|
[11123] | 762 | (see Chapters 13-16 of \cite{griffies_bk04}) for a discussion of neutral physics in $z$-models, |
---|
[10354] | 763 | as well as \autoref{sec:LDF_slp} in this document for treatment in \NEMO). |
---|
[2282] | 764 | |
---|
[13463] | 765 | The range over which \zstar\ varies is time independent $-H \leq \zstar \leq 0$. |
---|
[11335] | 766 | Hence, all cells remain nonvanishing, so long as the surface height maintains $\eta > -H$. |
---|
[13463] | 767 | This is a minor constraint relative to that encountered on the surface height when |
---|
| 768 | using $s = z$ or $s = z - \eta$. |
---|
[2282] | 769 | |
---|
[13463] | 770 | Because \zstar\ has a time independent range, all grid cells have static increments ds, |
---|
[11335] | 771 | and the sum of the vertical increments yields the time independent ocean depth. %k ds = H. |
---|
[13463] | 772 | The \zstar\ coordinate is therefore invisible to undulations of the free surface, |
---|
[10354] | 773 | since it moves along with the free surface. |
---|
[13463] | 774 | This property means that no spurious vertical transport is induced across |
---|
| 775 | surfaces of constant \zstar\ by the motion of external gravity waves. |
---|
[11335] | 776 | Such spurious transport can be a problem in z-models, especially those with tidal forcing. |
---|
[13463] | 777 | Quite generally, |
---|
| 778 | the time independent range for the \zstar\ coordinate is a very convenient property that |
---|
| 779 | allows for a nearly arbitrary vertical resolution even in |
---|
| 780 | the presence of large amplitude fluctuations of the surface height, again so long as $\eta > -H$. |
---|
[2282] | 781 | %end MOM doc %%% |
---|
| 782 | |
---|
[13463] | 783 | %% ================================================================================================= |
---|
[9393] | 784 | \subsection{Curvilinear terrain-following \textit{s}--coordinate} |
---|
[13463] | 785 | \label{subsec:MB_sco} |
---|
[707] | 786 | |
---|
[10354] | 787 | Several important aspects of the ocean circulation are influenced by bottom topography. |
---|
[13463] | 788 | Of course, the most important is that bottom topography determines deep ocean sub-basins, barriers, |
---|
| 789 | sills and channels that strongly constrain the path of water masses, but more subtle effects exist. |
---|
| 790 | For example, |
---|
| 791 | the topographic $\beta$-effect is usually larger than the planetary one along continental slopes. |
---|
[10354] | 792 | Topographic Rossby waves can be excited and can interact with the mean current. |
---|
[13463] | 793 | In the $z$-coordinate system presented in the previous section (\autoref{sec:MB_zco}), |
---|
[10501] | 794 | $z$-surfaces are geopotential surfaces. |
---|
[10354] | 795 | The bottom topography is discretised by steps. |
---|
| 796 | This often leads to a misrepresentation of a gradually sloping bottom and to |
---|
| 797 | large localized depth gradients associated with large localized vertical velocities. |
---|
| 798 | The response to such a velocity field often leads to numerical dispersion effects. |
---|
[13463] | 799 | One solution to strongly reduce this error is to use a partial step representation of |
---|
| 800 | bottom topography instead of a full step one \cite{pacanowski.gnanadesikan_MWR98}. |
---|
[10501] | 801 | Another solution is to introduce a terrain-following coordinate system (hereafter $s$-coordinate). |
---|
[707] | 802 | |
---|
[13463] | 803 | The $s$-coordinate avoids the discretisation error in the depth field since |
---|
| 804 | the layers of computation are gradually adjusted with depth to the ocean bottom. |
---|
| 805 | Relatively small topographic features as well as |
---|
| 806 | gentle, large-scale slopes of the sea floor in the deep ocean, |
---|
| 807 | which would be ignored in typical $z$-model applications with |
---|
| 808 | the largest grid spacing at greatest depths, |
---|
[10354] | 809 | can easily be represented (with relatively low vertical resolution). |
---|
[13463] | 810 | A terrain-following model (hereafter $s$-model) also facilitates |
---|
| 811 | the modelling of the boundary layer flows over a large depth range, |
---|
| 812 | which in the framework of the $z$-model would require high vertical resolution over |
---|
[10354] | 813 | the whole depth range. |
---|
[13463] | 814 | Moreover, |
---|
| 815 | with a $s$-coordinate it is possible, at least in principle, to have the bottom and the sea surface as |
---|
[10354] | 816 | the only boundaries of the domain (no more lateral boundary condition to specify). |
---|
[13463] | 817 | Nevertheless, a $s$-coordinate also has its drawbacks. |
---|
| 818 | Perfectly adapted to a homogeneous ocean, |
---|
[10354] | 819 | it has strong limitations as soon as stratification is introduced. |
---|
| 820 | The main two problems come from the truncation error in the horizontal pressure gradient and |
---|
| 821 | a possibly increased diapycnal diffusion. |
---|
[13463] | 822 | The horizontal pressure force in $s$-coordinate consists of two terms (see \autoref{apdx:SCOORD}), |
---|
[10414] | 823 | \begin{equation} |
---|
[13463] | 824 | \label{eq:MB_p_sco} |
---|
[11335] | 825 | \nabla p |_z = \nabla p |_s - \frac{1}{e_3} \pd[p]{s} \nabla z |_s |
---|
[707] | 826 | \end{equation} |
---|
| 827 | |
---|
[13463] | 828 | The second term in \autoref{eq:MB_p_sco} depends on the tilt of the coordinate surface and |
---|
[11335] | 829 | leads to a truncation error that is not present in a $z$-model. |
---|
[13463] | 830 | In the special case of a $\sigma$-coordinate |
---|
| 831 | (i.e. a depth-normalised coordinate system $\sigma = z/H$), |
---|
| 832 | \citet{haney_JPO91} and \citet{beckmann.haidvogel_JPO93} have given estimates of |
---|
| 833 | the magnitude of this truncation error. |
---|
| 834 | It depends on topographic slope, stratification, horizontal and vertical resolution, |
---|
| 835 | the equation of state, and the finite difference scheme. |
---|
[10354] | 836 | This error limits the possible topographic slopes that a model can handle at |
---|
| 837 | a given horizontal and vertical resolution. |
---|
| 838 | This is a severe restriction for large-scale applications using realistic bottom topography. |
---|
[13463] | 839 | The large-scale slopes require high horizontal resolution, |
---|
| 840 | and the computational cost becomes prohibitive. |
---|
[10354] | 841 | This problem can be at least partially overcome by mixing $s$-coordinate and |
---|
[13463] | 842 | step-like representation of bottom topography |
---|
| 843 | \citep{gerdes_JGR93*a,gerdes_JGR93*b,madec.delecluse.ea_JPO96}. |
---|
[10354] | 844 | However, the definition of the model domain vertical coordinate becomes then a non-trivial thing for |
---|
| 845 | a realistic bottom topography: |
---|
[13463] | 846 | an envelope topography is defined in $s$-coordinate on which |
---|
| 847 | a full or partial step bottom topography is then applied in order to |
---|
| 848 | adjust the model depth to the observed one (see \autoref{subsec:DOM_zgr}). |
---|
[707] | 849 | |
---|
[13463] | 850 | For numerical reasons a minimum of diffusion is required along |
---|
| 851 | the coordinate surfaces of any finite difference model. |
---|
[10354] | 852 | It causes spurious diapycnal mixing when coordinate surfaces do not coincide with isoneutral surfaces. |
---|
| 853 | This is the case for a $z$-model as well as for a $s$-model. |
---|
[13463] | 854 | However, density varies more strongly on $s$-surfaces than on horizontal surfaces in |
---|
| 855 | regions of large topographic slopes, |
---|
| 856 | implying larger diapycnal diffusion in a $s$-model than in a $z$-model. |
---|
| 857 | Whereas such a diapycnal diffusion in a $z$-model tends to |
---|
| 858 | weaken horizontal density (pressure) gradients and thus the horizontal circulation, |
---|
| 859 | it usually reinforces these gradients in a $s$-model, creating spurious circulation. |
---|
| 860 | For example, imagine an isolated bump of topography in |
---|
| 861 | an ocean at rest with a horizontally uniform stratification. |
---|
[10354] | 862 | Spurious diffusion along $s$-surfaces will induce a bump of isoneutral surfaces over the topography, |
---|
| 863 | and thus will generate there a baroclinic eddy. |
---|
| 864 | In contrast, the ocean will stay at rest in a $z$-model. |
---|
[13463] | 865 | As for the truncation error, the problem can be reduced by |
---|
| 866 | introducing the terrain-following coordinate below the strongly stratified portion of the water column |
---|
| 867 | (\ie\ the main thermocline) \citep{madec.delecluse.ea_JPO96}. |
---|
| 868 | An alternate solution consists of rotating the lateral diffusive tensor to |
---|
| 869 | geopotential or to isoneutral surfaces (see \autoref{subsec:MB_ldf}). |
---|
[10354] | 870 | Unfortunately, the slope of isoneutral surfaces relative to the $s$-surfaces can very large, |
---|
[13463] | 871 | strongly exceeding the stability limit of such a operator when it is discretized |
---|
| 872 | (see \autoref{chap:LDF}). |
---|
[707] | 873 | |
---|
[13463] | 874 | The $s$-coordinates introduced here \citep{lott.madec.ea_OM90,madec.delecluse.ea_JPO96} |
---|
| 875 | differ mainly in two aspects from similar models: |
---|
| 876 | it allows a representation of bottom topography with |
---|
| 877 | mixed full or partial step-like/terrain following topography; |
---|
[10354] | 878 | It also offers a completely general transformation, $s=s(i,j,z)$ for the vertical coordinate. |
---|
[707] | 879 | |
---|
[13463] | 880 | %% ================================================================================================= |
---|
| 881 | \subsection{Curvilinear \ztilde-coordinate} |
---|
| 882 | \label{subsec:MB_zco_tilde} |
---|
[707] | 883 | |
---|
[13463] | 884 | The \ztilde-coordinate has been developed by \citet{leclair.madec_OM11}. |
---|
| 885 | It is available in \NEMO\ since the version 3.4 and is more robust in version 4.0 than previously. |
---|
[10354] | 886 | Nevertheless, it is currently not robust enough to be used in all possible configurations. |
---|
| 887 | Its use is therefore not recommended. |
---|
[707] | 888 | |
---|
[13463] | 889 | %% ================================================================================================= |
---|
[9393] | 890 | \section{Subgrid scale physics} |
---|
[13463] | 891 | \label{sec:MB_zdf_ldf} |
---|
[707] | 892 | |
---|
[13463] | 893 | The hydrostatic primitive equations describe the behaviour of a geophysical fluid at |
---|
| 894 | space and time scales larger than a few kilometres in the horizontal, |
---|
| 895 | a few meters in the vertical and a few minutes. |
---|
| 896 | They are usually solved at larger scales: the specified grid spacing and time step of |
---|
| 897 | the numerical model. |
---|
[10501] | 898 | The effects of smaller scale motions (coming from the advective terms in the Navier-Stokes equations) |
---|
| 899 | must be represented entirely in terms of large-scale patterns to close the equations. |
---|
[10354] | 900 | These effects appear in the equations as the divergence of turbulent fluxes |
---|
[13463] | 901 | (\ie\ fluxes associated with the mean correlation of small scale perturbations). |
---|
[10354] | 902 | Assuming a turbulent closure hypothesis is equivalent to choose a formulation for these fluxes. |
---|
| 903 | It is usually called the subgrid scale physics. |
---|
| 904 | It must be emphasized that this is the weakest part of the primitive equations, |
---|
| 905 | but also one of the most important for long-term simulations as |
---|
| 906 | small scale processes \textit{in fine} balance the surface input of kinetic energy and heat. |
---|
[707] | 907 | |
---|
[13463] | 908 | The control exerted by gravity on the flow induces a strong anisotropy between |
---|
| 909 | the lateral and vertical motions. |
---|
[10501] | 910 | Therefore subgrid-scale physics \textbf{D}$^{\vect U}$, $D^{S}$ and $D^{T}$ in |
---|
[13463] | 911 | \autoref{eq:MB_PE_dyn}, \autoref{eq:MB_PE_tra_T} and \autoref{eq:MB_PE_tra_S} are divided into |
---|
| 912 | a lateral part \textbf{D}$^{l \vect U}$, $D^{l S}$ and $D^{l T}$ and |
---|
[10501] | 913 | a vertical part \textbf{D}$^{v \vect U}$, $D^{v S}$ and $D^{v T}$. |
---|
[13463] | 914 | The formulation of these terms and their underlying physics are briefly discussed in |
---|
| 915 | the next two subsections. |
---|
[707] | 916 | |
---|
[13463] | 917 | %% ================================================================================================= |
---|
[9393] | 918 | \subsection{Vertical subgrid scale physics} |
---|
[13463] | 919 | \label{subsec:MB_zdf} |
---|
[707] | 920 | |
---|
[13463] | 921 | The model resolution is always larger than the scale at which |
---|
| 922 | the major sources of vertical turbulence occur (shear instability, internal wave breaking...). |
---|
[10354] | 923 | Turbulent motions are thus never explicitly solved, even partially, but always parameterized. |
---|
[13463] | 924 | The vertical turbulent fluxes are assumed to depend linearly on |
---|
| 925 | the gradients of large-scale quantities (for example, |
---|
| 926 | the turbulent heat flux is given by $\overline{T' w'} = -A^{v T} \partial_z \overline T$, |
---|
[10501] | 927 | where $A^{v T}$ is an eddy coefficient). |
---|
[10354] | 928 | This formulation is analogous to that of molecular diffusion and dissipation. |
---|
| 929 | This is quite clearly a necessary compromise: considering only the molecular viscosity acting on |
---|
| 930 | large scale severely underestimates the role of turbulent diffusion and dissipation, |
---|
| 931 | while an accurate consideration of the details of turbulent motions is simply impractical. |
---|
| 932 | The resulting vertical momentum and tracer diffusive operators are of second order: |
---|
[10414] | 933 | \begin{equation} |
---|
[13463] | 934 | \label{eq:MB_zdf} |
---|
[10501] | 935 | \begin{gathered} |
---|
[13463] | 936 | \vect D^{v \vect U} = \pd[]{z} \lt( A^{vm} \pd[\vect U_h]{z} \rt) \text{,} \ |
---|
| 937 | D^{vT} = \pd[]{z} \lt( A^{vT} \pd[T]{z} \rt) \ \text{and} \ |
---|
[10501] | 938 | D^{vS} = \pd[]{z} \lt( A^{vT} \pd[S]{z} \rt) |
---|
| 939 | \end{gathered} |
---|
[707] | 940 | \end{equation} |
---|
[13463] | 941 | where $A^{vm}$ and $A^{vT}$ are the vertical eddy viscosity and diffusivity coefficients, |
---|
| 942 | respectively. |
---|
[10354] | 943 | At the sea surface and at the bottom, turbulent fluxes of momentum, heat and salt must be specified |
---|
| 944 | (see \autoref{chap:SBC} and \autoref{chap:ZDF} and \autoref{sec:TRA_bbl}). |
---|
| 945 | All the vertical physics is embedded in the specification of the eddy coefficients. |
---|
| 946 | They can be assumed to be either constant, or function of the local fluid properties |
---|
[13463] | 947 | (\eg\ Richardson number, Brunt-V\"{a}is\"{a}l\"{a} frequency, distance from the boundary ...), |
---|
[10354] | 948 | or computed from a turbulent closure model. |
---|
[13463] | 949 | The choices available in \NEMO\ are discussed in \autoref{chap:ZDF}). |
---|
[707] | 950 | |
---|
[13463] | 951 | %% ================================================================================================= |
---|
[9393] | 952 | \subsection{Formulation of the lateral diffusive and viscous operators} |
---|
[13463] | 953 | \label{subsec:MB_ldf} |
---|
[707] | 954 | |
---|
[10354] | 955 | Lateral turbulence can be roughly divided into a mesoscale turbulence associated with eddies |
---|
| 956 | (which can be solved explicitly if the resolution is sufficient since |
---|
| 957 | their underlying physics are included in the primitive equations), |
---|
[13463] | 958 | and a sub mesoscale turbulence which is never explicitly solved even partially, |
---|
| 959 | but always parameterized. |
---|
| 960 | The formulation of lateral eddy fluxes depends on whether |
---|
| 961 | the mesoscale is below or above the grid-spacing (\ie\ the model is eddy-resolving or not). |
---|
[707] | 962 | |
---|
[10354] | 963 | In non-eddy-resolving configurations, the closure is similar to that used for the vertical physics. |
---|
[13463] | 964 | The lateral turbulent fluxes are assumed to depend linearly on |
---|
| 965 | the lateral gradients of large-scale quantities. |
---|
[10354] | 966 | The resulting lateral diffusive and dissipative operators are of second order. |
---|
[13463] | 967 | Observations show that |
---|
| 968 | lateral mixing induced by mesoscale turbulence tends to be along isopycnal surfaces |
---|
[11123] | 969 | (or more precisely neutral surfaces \cite{mcdougall_JPO87}) rather than across them. |
---|
[13463] | 970 | As the slope of neutral surfaces is small in the ocean, |
---|
| 971 | a common approximation is to assume that the ``lateral'' direction is the horizontal, |
---|
| 972 | \ie\ the lateral mixing is performed along geopotential surfaces. |
---|
[10354] | 973 | This leads to a geopotential second order operator for lateral subgrid scale physics. |
---|
[13463] | 974 | This assumption can be relaxed: |
---|
| 975 | the eddy-induced turbulent fluxes can be better approached by assuming that |
---|
[10354] | 976 | they depend linearly on the gradients of large-scale quantities computed along neutral surfaces. |
---|
| 977 | In such a case, the diffusive operator is an isoneutral second order operator and |
---|
| 978 | it has components in the three space directions. |
---|
[13463] | 979 | However, both horizontal and isoneutral operators have no effect on |
---|
| 980 | mean (\ie\ large scale) potential energy whereas |
---|
[10354] | 981 | potential energy is a main source of turbulence (through baroclinic instabilities). |
---|
[11335] | 982 | \citet{gent.mcwilliams_JPO90} proposed a parameterisation of mesoscale eddy-induced turbulence which |
---|
[10354] | 983 | associates an eddy-induced velocity to the isoneutral diffusion. |
---|
| 984 | Its mean effect is to reduce the mean potential energy of the ocean. |
---|
[13463] | 985 | This leads to a formulation of lateral subgrid-scale physics made up of |
---|
| 986 | an isoneutral second order operator and an eddy induced advective part. |
---|
[10354] | 987 | In all these lateral diffusive formulations, |
---|
| 988 | the specification of the lateral eddy coefficients remains the problematic point as |
---|
[13463] | 989 | there is no really satisfactory formulation of these coefficients as |
---|
| 990 | a function of large-scale features. |
---|
[707] | 991 | |
---|
[10354] | 992 | In eddy-resolving configurations, a second order operator can be used, |
---|
| 993 | but usually the more scale selective biharmonic operator is preferred as |
---|
| 994 | the grid-spacing is usually not small enough compared to the scale of the eddies. |
---|
| 995 | The role devoted to the subgrid-scale physics is to dissipate the energy that |
---|
| 996 | cascades toward the grid scale and thus to ensure the stability of the model while |
---|
| 997 | not interfering with the resolved mesoscale activity. |
---|
| 998 | Another approach is becoming more and more popular: |
---|
[13463] | 999 | instead of specifying explicitly a sub-grid scale term in |
---|
| 1000 | the momentum and tracer time evolution equations, |
---|
[11335] | 1001 | one uses an advective scheme which is diffusive enough to maintain the model stability. |
---|
[13463] | 1002 | It must be emphasised that then, |
---|
| 1003 | all the sub-grid scale physics is included in the formulation of the advection scheme. |
---|
[707] | 1004 | |
---|
[10354] | 1005 | All these parameterisations of subgrid scale physics have advantages and drawbacks. |
---|
[13463] | 1006 | They are not all available in \NEMO. |
---|
| 1007 | For active tracers (temperature and salinity) the main ones are: |
---|
[10354] | 1008 | Laplacian and bilaplacian operators acting along geopotential or iso-neutral surfaces, |
---|
[11123] | 1009 | \citet{gent.mcwilliams_JPO90} parameterisation, and various slightly diffusive advection schemes. |
---|
[13463] | 1010 | For momentum, the main ones are: |
---|
| 1011 | Laplacian and bilaplacian operators acting along geopotential surfaces, |
---|
[10354] | 1012 | and UBS advection schemes when flux form is chosen for the momentum advection. |
---|
[707] | 1013 | |
---|
[13463] | 1014 | %% ================================================================================================= |
---|
[9393] | 1015 | \subsubsection{Lateral laplacian tracer diffusive operator} |
---|
[707] | 1016 | |
---|
[13463] | 1017 | The lateral Laplacian tracer diffusive operator is defined by (see \autoref{apdx:DIFFOPERS}): |
---|
[10414] | 1018 | \begin{equation} |
---|
[13463] | 1019 | \label{eq:MB_iso_tensor} |
---|
| 1020 | D^{lT} = \nabla \vect . \lt( A^{lT} \; \Re \; \nabla T \rt) \quad \text{with} \quad \Re = |
---|
| 1021 | \begin{pmatrix} |
---|
| 1022 | 1 & 0 & -r_1 \\ |
---|
| 1023 | 0 & 1 & -r_2 \\ |
---|
| 1024 | -r_1 & -r_2 & r_1^2 + r_2^2 \\ |
---|
| 1025 | \end{pmatrix} |
---|
[707] | 1026 | \end{equation} |
---|
[10501] | 1027 | where $r_1$ and $r_2$ are the slopes between the surface along which the diffusive operator acts and |
---|
[13463] | 1028 | the model level (\eg\ $z$- or $s$-surfaces). |
---|
| 1029 | Note that the formulation \autoref{eq:MB_iso_tensor} is exact for |
---|
[10354] | 1030 | the rotation between geopotential and $s$-surfaces, |
---|
| 1031 | while it is only an approximation for the rotation between isoneutral and $z$- or $s$-surfaces. |
---|
[13463] | 1032 | Indeed, in the latter case, two assumptions are made to simplify \autoref{eq:MB_iso_tensor} |
---|
| 1033 | \citep{cox_OM87}. |
---|
| 1034 | First, the horizontal contribution of the dianeutral mixing is neglected since |
---|
| 1035 | the ratio between iso and dia-neutral diffusive coefficients is known to be |
---|
| 1036 | several orders of magnitude smaller than unity. |
---|
[10354] | 1037 | Second, the two isoneutral directions of diffusion are assumed to be independent since |
---|
[13463] | 1038 | the slopes are generally less than $10^{-2}$ in the ocean (see \autoref{apdx:DIFFOPERS}). |
---|
[707] | 1039 | |
---|
[10354] | 1040 | For \textit{iso-level} diffusion, $r_1$ and $r_2 $ are zero. |
---|
[10501] | 1041 | $\Re$ reduces to the identity in the horizontal direction, no rotation is applied. |
---|
[6140] | 1042 | |
---|
[10354] | 1043 | For \textit{geopotential} diffusion, |
---|
| 1044 | $r_1$ and $r_2 $ are the slopes between the geopotential and computational surfaces: |
---|
[13463] | 1045 | they are equal to $\sigma_1$ and $\sigma_2$, respectively (see \autoref{eq:MB_sco_slope}). |
---|
[707] | 1046 | |
---|
[13463] | 1047 | For \textit{isoneutral} diffusion $r_1$ and $r_2$ are the slopes between |
---|
| 1048 | the isoneutral and computational surfaces. |
---|
[10354] | 1049 | Therefore, they are different quantities, but have similar expressions in $z$- and $s$-coordinates. |
---|
| 1050 | In $z$-coordinates: |
---|
[10414] | 1051 | \begin{equation} |
---|
[13463] | 1052 | \label{eq:MB_iso_slopes} |
---|
| 1053 | r_1 = \frac{e_3}{e_1} \lt( \pd[\rho]{i} \rt) \lt( \pd[\rho]{k} \rt)^{-1} \quad |
---|
| 1054 | r_2 = \frac{e_3}{e_2} \lt( \pd[\rho]{j} \rt) \lt( \pd[\rho]{k} \rt)^{-1} |
---|
[707] | 1055 | \end{equation} |
---|
[6289] | 1056 | while in $s$-coordinates $\pd[]{k}$ is replaced by $\pd[]{s}$. |
---|
[707] | 1057 | |
---|
[13463] | 1058 | %% ================================================================================================= |
---|
[3294] | 1059 | \subsubsection{Eddy induced velocity} |
---|
[10501] | 1060 | |
---|
[11123] | 1061 | When the \textit{eddy induced velocity} parametrisation (eiv) \citep{gent.mcwilliams_JPO90} is used, |
---|
[1224] | 1062 | an additional tracer advection is introduced in combination with the isoneutral diffusion of tracers: |
---|
[10414] | 1063 | \[ |
---|
[13463] | 1064 | % \label{eq:MB_iso+eiv} |
---|
[10501] | 1065 | D^{lT} = \nabla \cdot \lt( A^{lT} \; \Re \; \nabla T \rt) + \nabla \cdot \lt( \vect U^\ast \, T \rt) |
---|
[10414] | 1066 | \] |
---|
[10501] | 1067 | where $ \vect U^\ast = \lt( u^\ast,v^\ast,w^\ast \rt)$ is a non-divergent, |
---|
[1224] | 1068 | eddy-induced transport velocity. This velocity field is defined by: |
---|
[13463] | 1069 | \[ |
---|
| 1070 | % \label{eq:MB_eiv} |
---|
| 1071 | u^\ast = \frac{1}{e_3} \pd[]{k} \lt( A^{eiv} \; \tilde{r}_1 \rt) \quad |
---|
| 1072 | v^\ast = \frac{1}{e_3} \pd[]{k} \lt( A^{eiv} \; \tilde{r}_2 \rt) \quad |
---|
[10544] | 1073 | w^\ast = - \frac{1}{e_1 e_2} \lt[ \pd[]{i} \lt( A^{eiv} \; e_2 \, \tilde{r}_1 \rt) |
---|
[10501] | 1074 | + \pd[]{j} \lt( A^{eiv} \; e_1 \, \tilde{r}_2 \rt) \rt] |
---|
[13463] | 1075 | \] |
---|
[10354] | 1076 | where $A^{eiv}$ is the eddy induced velocity coefficient |
---|
| 1077 | (or equivalently the isoneutral thickness diffusivity coefficient), |
---|
[13463] | 1078 | and $\tilde r_1$ and $\tilde r_2$ are the slopes between |
---|
| 1079 | isoneutral and \textit{geopotential} surfaces. |
---|
| 1080 | Their values are thus independent of the vertical coordinate, |
---|
| 1081 | but their expression depends on the coordinate: |
---|
| 1082 | \begin{equation} |
---|
| 1083 | \label{eq:MB_slopes_eiv} |
---|
[10414] | 1084 | \tilde{r}_n = |
---|
[10501] | 1085 | \begin{cases} |
---|
| 1086 | r_n & \text{in $z$-coordinate} \\ |
---|
| 1087 | r_n + \sigma_n & \text{in \zstar- and $s$-coordinates} |
---|
| 1088 | \end{cases} |
---|
| 1089 | \quad \text{where~} n = 1, 2 |
---|
[13463] | 1090 | \end{equation} |
---|
[707] | 1091 | |
---|
[10354] | 1092 | The normal component of the eddy induced velocity is zero at all the boundaries. |
---|
[13463] | 1093 | This can be achieved in a model by tapering either the eddy coefficient or the slopes to zero in |
---|
| 1094 | the vicinity of the boundaries. |
---|
| 1095 | The latter strategy is used in \NEMO\ (cf. \autoref{chap:LDF}). |
---|
[707] | 1096 | |
---|
[13463] | 1097 | %% ================================================================================================= |
---|
[6289] | 1098 | \subsubsection{Lateral bilaplacian tracer diffusive operator} |
---|
[707] | 1099 | |
---|
[6289] | 1100 | The lateral bilaplacian tracer diffusive operator is defined by: |
---|
[10414] | 1101 | \[ |
---|
[13463] | 1102 | % \label{eq:MB_bilapT} |
---|
[10501] | 1103 | D^{lT}= - \Delta \; (\Delta T) \quad \text{where} \quad |
---|
| 1104 | \Delta \bullet = \nabla \lt( \sqrt{B^{lT}} \; \Re \; \nabla \bullet \rt) |
---|
[10414] | 1105 | \] |
---|
[13463] | 1106 | It is the Laplacian operator given by \autoref{eq:MB_iso_tensor} applied twice with |
---|
[10501] | 1107 | the harmonic eddy diffusion coefficient set to the square root of the biharmonic one. |
---|
[707] | 1108 | |
---|
[13463] | 1109 | %% ================================================================================================= |
---|
[6289] | 1110 | \subsubsection{Lateral Laplacian momentum diffusive operator} |
---|
[707] | 1111 | |
---|
[10354] | 1112 | The Laplacian momentum diffusive operator along $z$- or $s$-surfaces is found by |
---|
[13463] | 1113 | applying \autoref{eq:MB_lap_vector} to the horizontal velocity vector (see \autoref{apdx:DIFFOPERS}): |
---|
[10501] | 1114 | \begin{align*} |
---|
[13463] | 1115 | % \label{eq:MB_lapU} |
---|
[10501] | 1116 | \vect D^{l \vect U} &= \nabla_h \big( A^{lm} \chi \big) |
---|
| 1117 | - \nabla_h \times \big( A^{lm} \, \zeta \; \vect k \big) \\ |
---|
| 1118 | &= \lt( \frac{1}{e_1} \pd[ \lt( A^{lm} \chi \rt) ]{i} \rt. |
---|
[13463] | 1119 | - \frac{1}{e_2 e_3} \pd[ \lt( A^{lm} \; e_3 \zeta \rt) ]{j} , |
---|
[10501] | 1120 | \frac{1}{e_2} \pd[ \lt( A^{lm} \chi \rt) ]{j} |
---|
| 1121 | \lt. + \frac{1}{e_1 e_3} \pd[ \lt( A^{lm} \; e_3 \zeta \rt) ]{i} \rt) |
---|
| 1122 | \end{align*} |
---|
[707] | 1123 | |
---|
[13463] | 1124 | Such a formulation ensures a complete separation between |
---|
| 1125 | the vorticity and horizontal divergence fields (see \autoref{apdx:INVARIANTS}). |
---|
[10354] | 1126 | Unfortunately, it is only available in \textit{iso-level} direction. |
---|
| 1127 | When a rotation is required |
---|
[13463] | 1128 | (\ie\ geopotential diffusion in $s$-coordinates or |
---|
| 1129 | isoneutral diffusion in both $z$- and $s$-coordinates), |
---|
| 1130 | the $u$ and $v$-fields are considered as independent scalar fields, |
---|
| 1131 | so that the diffusive operator is given by: |
---|
| 1132 | \[ |
---|
| 1133 | % \label{eq:MB_lapU_iso} |
---|
| 1134 | D_u^{l \vect U} = \nabla . \lt( A^{lm} \; \Re \; \nabla u \rt) \quad |
---|
[10501] | 1135 | D_v^{l \vect U} = \nabla . \lt( A^{lm} \; \Re \; \nabla v \rt) |
---|
[13463] | 1136 | \] |
---|
| 1137 | where $\Re$ is given by \autoref{eq:MB_iso_tensor}. |
---|
[10354] | 1138 | It is the same expression as those used for diffusive operator on tracers. |
---|
| 1139 | It must be emphasised that such a formulation is only exact in a Cartesian coordinate system, |
---|
[13463] | 1140 | \ie\ on a $f$- or $\beta$-plane, not on the sphere. |
---|
[10354] | 1141 | It is also a very good approximation in vicinity of the Equator in |
---|
[11123] | 1142 | a geographical coordinate system \citep{lengaigne.madec.ea_JGR03}. |
---|
[707] | 1143 | |
---|
[13463] | 1144 | %% ================================================================================================= |
---|
[11335] | 1145 | \subsubsection{Lateral bilaplacian momentum diffusive operator} |
---|
[707] | 1146 | |
---|
[13463] | 1147 | As for tracers, |
---|
| 1148 | the bilaplacian order momentum diffusive operator is a re-entering Laplacian operator with |
---|
[10354] | 1149 | the harmonic eddy diffusion coefficient set to the square root of the biharmonic one. |
---|
| 1150 | Nevertheless it is currently not available in the iso-neutral case. |
---|
[707] | 1151 | |
---|
[13463] | 1152 | \subinc{\input{../../global/epilogue}} |
---|
[10414] | 1153 | |
---|
[6997] | 1154 | \end{document} |
---|