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