[707] | 1 | % ================================================================ |
---|
| 2 | % Chapter 1 Ñ Model Basics |
---|
| 3 | % ================================================================ |
---|
| 4 | |
---|
| 5 | \chapter{Model basics} |
---|
| 6 | \label{PE} |
---|
| 7 | \minitoc |
---|
| 8 | |
---|
| 9 | |
---|
| 10 | % ================================================================ |
---|
| 11 | % Primitive Equations |
---|
| 12 | % ================================================================ |
---|
| 13 | \section{Primitive Equations} |
---|
| 14 | \label{PE_PE} |
---|
| 15 | |
---|
| 16 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 17 | % Vector Invariant Formulation |
---|
| 18 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 19 | |
---|
| 20 | \subsection{Vector Invariant Formulation} |
---|
| 21 | \label{PE_Vector} |
---|
| 22 | |
---|
| 23 | |
---|
[996] | 24 | The ocean is a fluid that can be described to a good approximation by the primitive |
---|
| 25 | equations, $i.e.$ the Navier-Stokes equations along with a nonlinear equation of |
---|
| 26 | state which couples the two active tracers (temperature and salinity) to the fluid |
---|
| 27 | velocity, plus the following additional assumptions made from scale considerations: |
---|
[707] | 28 | |
---|
[1224] | 29 | \textit{(1) spherical earth approximation: }the geopotential surfaces are assumed to |
---|
| 30 | be spheres so that gravity (local vertical) is parallel to the earth's radius |
---|
[707] | 31 | |
---|
| 32 | \textit{(2) thin-shell approximation: }the ocean depth is neglected compared to the earth's radius |
---|
| 33 | |
---|
[1224] | 34 | \textit{(3) turbulent closure hypothesis: }the turbulent fluxes (which represent the effect |
---|
| 35 | of small scale processes on the large-scale) are expressed in terms of large-scale features |
---|
[707] | 36 | |
---|
[1224] | 37 | \textit{(4) Boussinesq hypothesis:} density variations are neglected except in their |
---|
| 38 | contribution to the buoyancy force |
---|
[707] | 39 | |
---|
[1224] | 40 | \textit{(5) Hydrostatic hypothesis: }the vertical momentum equation is reduced to a |
---|
| 41 | balance between the vertical pressure gradient and the buoyancy force (this removes |
---|
| 42 | convective processes from the initial Navier-Stokes equations and so convective processes |
---|
| 43 | must be parameterized instead) |
---|
[707] | 44 | |
---|
[1224] | 45 | \textit{(6) Incompressibility hypothesis: }the three dimensional divergence of the velocity |
---|
| 46 | vector is assumed to be zero. |
---|
[707] | 47 | |
---|
[1224] | 48 | Because the gravitational force is so dominant in the equations of large-scale motions, |
---|
| 49 | it is useful to choose an orthogonal set of unit vectors (\textbf{i},\textbf{j},\textbf{k}) linked |
---|
| 50 | to the earth such that \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are two |
---|
| 51 | vectors orthogonal to \textbf{k}, $i.e.$ tangent to the geopotential surfaces. Let us define |
---|
| 52 | the following variables: \textbf{U} the vector velocity, $\textbf{U}=\textbf{U}_h + w\, \textbf{k}$ |
---|
| 53 | (the subscript $h$ denotes the local horizontal vector, $i.e.$ over the (\textbf{i},\textbf{j}) plane), |
---|
| 54 | $T$ the potential temperature, $S$ the salinity, \textit{$\rho $} the \textit{in situ} density. |
---|
| 55 | The vector invariant form of the primitive equations in the (\textbf{i},\textbf{j},\textbf{k}) |
---|
| 56 | vector system provides the following six equations (namely the momentum balance, the |
---|
| 57 | hydrostatic equilibrium, the incompressibility equation, the heat and salt conservation |
---|
| 58 | equations and an equation of state): |
---|
[707] | 59 | \begin{subequations} \label{Eq_PE} |
---|
| 60 | \begin{equation} \label{Eq_PE_dyn} |
---|
| 61 | \frac{\partial {\rm {\bf U}}_h }{\partial t}= |
---|
| 62 | -\left[ {\left( {\nabla \times {\rm {\bf U}}} \right)\times {\rm {\bf U}} |
---|
| 63 | +\frac{1}{2}\nabla \left( {{\rm {\bf U}}^2} \right)} \right]_h |
---|
| 64 | -f\;{\rm {\bf k}}\times {\rm {\bf U}}_h |
---|
[817] | 65 | -\frac{1}{\rho _o }\nabla _h p + {\rm {\bf D}}^{\rm {\bf U}} + {\rm {\bf F}}^{\rm {\bf U}} |
---|
[707] | 66 | \end{equation} |
---|
| 67 | \begin{equation} \label{Eq_PE_hydrostatic} |
---|
| 68 | \frac{\partial p }{\partial z} = - \rho \ g |
---|
| 69 | \end{equation} |
---|
| 70 | \begin{equation} \label{Eq_PE_continuity} |
---|
| 71 | \nabla \cdot {\bf U}= 0 |
---|
| 72 | \end{equation} |
---|
| 73 | \begin{equation} \label{Eq_PE_tra_T} |
---|
[817] | 74 | \frac{\partial T}{\partial t} = - \nabla \cdot \left( T \ \rm{\bf U} \right) + D^T + F^T |
---|
[707] | 75 | \end{equation} |
---|
| 76 | \begin{equation} \label{Eq_PE_tra_S} |
---|
[817] | 77 | \frac{\partial S}{\partial t} = - \nabla \cdot \left( S \ \rm{\bf U} \right) + D^S + F^S |
---|
[707] | 78 | \end{equation} |
---|
| 79 | \begin{equation} \label{Eq_PE_eos} |
---|
| 80 | \rho = \rho \left( T,S,p \right) |
---|
| 81 | \end{equation} |
---|
| 82 | \end{subequations} |
---|
[1224] | 83 | where $\nabla$ is the generalised derivative vector operator in $(\bf i,\bf j, \bf k)$ directions, |
---|
| 84 | $t$ is the time, $z$ is the vertical coordinate, $\rho $ is the \textit{in situ} density given by |
---|
| 85 | the equation of state (\ref{Eq_PE_eos}), $\rho_o$ is a reference density, $p$ the pressure, |
---|
| 86 | $f=2 \bf \Omega \cdot \bf k$ is the Coriolis acceleration (where $\bf \Omega$ is the Earth's |
---|
| 87 | angular velocity vector), and $g$ is the gravitational acceleration. |
---|
| 88 | ${\rm {\bf D}}^{\rm {\bf U}}$, $D^T$ and $D^S$ are the parameterisations of small-scale |
---|
| 89 | physics for momentum, temperature and salinity, and ${\rm {\bf F}}^{\rm {\bf U}}$, $F^T$ |
---|
| 90 | and $F^S$ surface forcing terms. Their nature and formulation are discussed in |
---|
| 91 | \S\ref{PE_zdf_ldf} and page \S\ref{PE_boundary_condition}. |
---|
[707] | 92 | |
---|
| 93 | . |
---|
| 94 | |
---|
| 95 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 96 | % Boundary condition |
---|
| 97 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 98 | \subsection{Boundary Conditions} |
---|
| 99 | \label{PE_boundary_condition} |
---|
| 100 | |
---|
[1224] | 101 | An ocean is bounded by complex coastlines, bottom topography at its base and an air-sea |
---|
| 102 | or ice-sea interface at its top. These boundaries can be defined by two surfaces, $z=-H(i,j)$ |
---|
| 103 | and $z=\eta(i,j,k,t)$, where $H$ is the depth of the ocean bottom and $\eta$ is the height |
---|
| 104 | of the sea surface. Both $H$ and $\eta$ are usually referenced to a given surface, $z=0$, |
---|
| 105 | chosen as a mean sea surface (Fig.~\ref{Fig_ocean_bc}). Through these two boundaries, |
---|
| 106 | the ocean can exchange fluxes of heat, fresh water, salt, and momentum with the solid earth, |
---|
| 107 | the continental margins, the sea ice and the atmosphere. However, some of these fluxes are |
---|
| 108 | so weak that even on climatic time scales of thousands of years they can be neglected. |
---|
| 109 | In the following, we briefly review the fluxes exchanged at the interfaces between the ocean |
---|
| 110 | and the other components of the earth system. |
---|
[707] | 111 | |
---|
| 112 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 113 | \begin{figure}[!ht] \label{Fig_ocean_bc} \begin{center} |
---|
[998] | 114 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_I_ocean_bc.pdf} |
---|
[1224] | 115 | \caption{The ocean is bounded by two surfaces, $z=-H(i,j)$ and $z=\eta(i,j,k,t)$, where $H$ |
---|
| 116 | is the depth of the sea floor and $\eta$ the height of the sea surface. Both $H$ and $\eta $ |
---|
| 117 | are referenced to $z=0$.} |
---|
[707] | 118 | \end{center} \end{figure} |
---|
| 119 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 120 | |
---|
[817] | 121 | |
---|
[707] | 122 | \begin{description} |
---|
[1224] | 123 | \item[Land - ocean interface:] the major flux between continental margins and the ocean is |
---|
| 124 | a mass exchange of fresh water through river runoff. Such an exchange modifies the sea |
---|
| 125 | surface salinity especially in the vicinity of major river mouths. It can be neglected for short |
---|
| 126 | range integrations but has to be taken into account for long term integrations as it influences |
---|
| 127 | the characteristics of water masses formed (especially at high latitudes). It is required in order |
---|
| 128 | to close the water cycle of the climate system. It is usually specified as a fresh water flux at |
---|
| 129 | the air-sea interface in the vicinity of river mouths. |
---|
| 130 | \item[Solid earth - ocean interface:] heat and salt fluxes through the sea floor are small, |
---|
| 131 | except in special areas of little extent. They are usually neglected in the model |
---|
| 132 | \footnote{In fact, it has been shown that the heat flux associated with the solid Earth cooling |
---|
| 133 | ($i.e.$the geothermal heating) is not negligible for the thermohaline circulation of the world |
---|
| 134 | ocean (see \ref{TRA_bbc}).}. |
---|
| 135 | The boundary condition is thus set to no flux of heat and salt across solid boundaries. |
---|
| 136 | For momentum, the situation is different. There is no flow across solid boundaries, |
---|
| 137 | $i.e.$ the velocity normal to the ocean bottom and coastlines is zero (in other words, |
---|
| 138 | the bottom velocity is parallel to solid boundaries). This kinematic boundary condition |
---|
| 139 | can be expressed as: |
---|
[707] | 140 | \begin{equation} \label{Eq_PE_w_bbc} |
---|
| 141 | w = -{\rm {\bf U}}_h \cdot \nabla _h \left( H \right) |
---|
| 142 | \end{equation} |
---|
[1224] | 143 | In addition, the ocean exchanges momentum with the earth through frictional processes. |
---|
| 144 | Such momentum transfer occurs at small scales in a boundary layer. It must be parameterized |
---|
| 145 | in terms of turbulent fluxes using bottom and/or lateral boundary conditions. Its specification |
---|
| 146 | depends on the nature of the physical parameterisation used for ${\rm {\bf D}}^{\rm {\bf U}}$ |
---|
| 147 | in \eqref{Eq_PE_dyn}. It is discussed in \S\ref{PE_zdf}, page~\pageref{PE_zdf}.% and Chap. III.6 to 9. |
---|
| 148 | \item[Atmosphere - ocean interface:] the kinematic surface condition plus the mass flux |
---|
| 149 | of fresh water PE (the precipitation minus evaporation budget) leads to: |
---|
[707] | 150 | \begin{equation} \label{Eq_PE_w_sbc} |
---|
| 151 | w = \frac{\partial \eta }{\partial t} |
---|
| 152 | + \left. {{\rm {\bf U}}_h } \right|_{z=\eta } \cdot \nabla _h \left( \eta \right) |
---|
| 153 | + P-E |
---|
| 154 | \end{equation} |
---|
[1224] | 155 | The dynamic boundary condition, neglecting the surface tension (which removes capillary |
---|
| 156 | waves from the system) leads to the continuity of pressure across the interface $z=\eta$. |
---|
| 157 | The atmosphere and ocean also exchange horizontal momentum (wind stress), and heat. |
---|
| 158 | \item[Sea ice - ocean interface:] the ocean and sea ice exchange heat, salt, fresh water |
---|
| 159 | and momentum. The sea surface temperature is constrained to be at the freezing point |
---|
| 160 | at the interface. Sea ice salinity is very low ($\sim4-6 \,psu$) compared to those of the |
---|
| 161 | ocean ($\sim34 \,psu$). The cycle of freezing/melting is associated with fresh water and |
---|
| 162 | salt fluxes that cannot be neglected. |
---|
[707] | 163 | \end{description} |
---|
| 164 | |
---|
| 165 | |
---|
| 166 | % ================================================================ |
---|
| 167 | % The Horizontal Pressure Gradient |
---|
| 168 | % ================================================================ |
---|
| 169 | \section{The Horizontal Pressure Gradient } |
---|
| 170 | \label{PE_hor_pg} |
---|
| 171 | |
---|
| 172 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 173 | % Pressure Formulation |
---|
| 174 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 175 | \subsection{Pressure Formulation} |
---|
| 176 | \label{PE_p_formulation} |
---|
| 177 | |
---|
[1224] | 178 | The total pressure at a given depth $z$ is composed of a surface pressure $p_s$ at a |
---|
| 179 | reference geopotential surface ($z=0$) and a hydrostatic pressure $p_h$ such that: |
---|
| 180 | $p(i,j,k,t)=p_s(i,j,t)+p_h(i,j,k,t)$. The latter is computed by integrating (\ref{Eq_PE_hydrostatic}), |
---|
| 181 | assuming that pressure in decibars can be approximated by depth in meters in (\ref{Eq_PE_eos}). |
---|
| 182 | The hydrostatic pressure is then given by: |
---|
[707] | 183 | \begin{equation} \label{Eq_PE_pressure} |
---|
| 184 | p_h \left( {i,j,z,t} \right) |
---|
[1224] | 185 | = \int_{\varsigma =z}^{\varsigma =0} {g\;\rho \left( {T,S,\varsigma} \right)\;d\varsigma } |
---|
[707] | 186 | \end{equation} |
---|
[1224] | 187 | Two strategies can be considered for the surface pressure term: $(a)$ introduce of a |
---|
| 188 | new variable $\eta$, the free-surface elevation, for which a prognostic equation can be |
---|
| 189 | established and solved; $(b)$ assume that the ocean surface is a rigid lid, on which the |
---|
| 190 | pressure (or its horizontal gradient) can be diagnosed. When the former strategy is used, |
---|
| 191 | one solution of the free-surface elevation consists of the excitation of external gravity waves. |
---|
| 192 | The flow is barotropic and the surface moves up and down with gravity as the restoring force. |
---|
| 193 | The phase speed of such waves is high (some hundreds of metres per second) so that |
---|
| 194 | the time step would have to be very short if they were present in the model. The latter |
---|
| 195 | strategy filters out these waves since the rigid lid approximation implies $\eta=0$, $i.e.$ |
---|
| 196 | the sea surface is the surface $z=0$. This well known approximation increases the surface |
---|
| 197 | wave speed to infinity and modifies certain other longwave dynamics ($e.g.$ barotropic |
---|
| 198 | Rossby or planetary waves). In the present release of \NEMO, both strategies are still available. |
---|
| 199 | They are further described in the next two sub-sections. |
---|
[707] | 200 | |
---|
| 201 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 202 | % Free Surface Formulation |
---|
| 203 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 204 | \subsection{Free Surface Formulation} |
---|
| 205 | \label{PE_free_surface} |
---|
| 206 | |
---|
[1224] | 207 | In the free surface formulation, a variable $\eta$, the sea-surface height, is introduced |
---|
| 208 | which describes the shape of the air-sea interface. This variable is solution of a |
---|
| 209 | prognostic equation which is established by forming the vertical average of the kinematic |
---|
| 210 | surface condition (\ref{Eq_PE_w_bbc}): |
---|
[707] | 211 | \begin{equation} \label{Eq_PE_ssh} |
---|
| 212 | \frac{\partial \eta }{\partial t}=-D+P-E |
---|
| 213 | \quad \text{where} \ |
---|
| 214 | D=\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\rm{\bf \overline{U}}}_h \,} \right] |
---|
| 215 | \end{equation} |
---|
| 216 | and using (\ref{Eq_PE_hydrostatic}) the surface pressure is given by: $p_s = \rho \, g \, \eta$. |
---|
| 217 | |
---|
[1224] | 218 | Allowing the air-sea interface to move introduces the external gravity waves (EGWs) |
---|
| 219 | as a class of solution of the primitive equations. These waves are barotropic because |
---|
| 220 | of hydrostatic assumption, and their phase speed is quite high. Their time scale is |
---|
| 221 | short with respect to the other processes described by the primitive equations. |
---|
[707] | 222 | |
---|
[1224] | 223 | Three choices can be made regarding the implementation of the free surface in the model, |
---|
| 224 | depending on the physical processes of interest. |
---|
[707] | 225 | |
---|
| 226 | $\bullet$ If one is interested in EGWs, in particular the tides and their interaction |
---|
[1224] | 227 | with the baroclinic structure of the ocean (internal waves) possibly in shallow seas, |
---|
| 228 | then a non linear free surface is the most appropriate. This means that no |
---|
| 229 | approximation is made in (\ref{Eq_PE_ssh}) and that the variation of the ocean |
---|
| 230 | volume is fully taken into account. Note that in order to study the fast time scales |
---|
| 231 | associated with EGWs it is necessary to minimize time filtering effects (use an |
---|
| 232 | explicit time scheme with very small time step, or a split-explicit scheme with |
---|
| 233 | reasonably small time step, see \S\ref{DYN_spg_exp} or \S\ref{DYN_spg_ts}. |
---|
[707] | 234 | |
---|
| 235 | $\bullet$ If one is not interested in EGW but rather sees them as high frequency |
---|
[817] | 236 | noise, it is possible to apply an explicit filter to slow down the fastest waves while |
---|
| 237 | not altering the slow barotropic Rossby waves. If further, an approximative conservation |
---|
| 238 | of heat and salt contents is sufficient for the problem solved, then it is |
---|
| 239 | sufficient to solve a linearized version of (\ref{Eq_PE_ssh}), which still allows |
---|
[707] | 240 | to take into account freshwater fluxes applied at the ocean surface \citep{Roullet2000}. |
---|
| 241 | |
---|
| 242 | $\bullet$ For process studies not involving external waves nor surface freshwater |
---|
| 243 | fluxes, it is possible to use the rigid lid approximation see (next |
---|
[817] | 244 | section). The ocean surface is then considered as a fixed surface, so that all |
---|
[707] | 245 | external waves are removed from the system. |
---|
| 246 | |
---|
[1224] | 247 | The filtering of EGWs in models with a free surface is usually a matter of discretisation |
---|
| 248 | of the temporal derivatives, using the time splitting method \citep{Killworth1991, Zhang1992} |
---|
| 249 | or the implicit scheme \citep{Dukowicz1994}. In \NEMO, we use a slightly different approach |
---|
| 250 | developed by \citet{Roullet2000}: the damping of EGWs is ensured by introducing an |
---|
| 251 | additional force in the momentum equation. \eqref{Eq_PE_dyn} becomes: |
---|
[707] | 252 | \begin{equation} \label{Eq_PE_flt} |
---|
| 253 | \frac{\partial {\rm {\bf U}}_h }{\partial t}= {\rm {\bf M}} |
---|
| 254 | - g \nabla \left( \tilde{\rho} \ \eta \right) |
---|
| 255 | - g \ T_c \nabla \left( \widetilde{\rho} \ \partial_t \eta \right) |
---|
| 256 | \end{equation} |
---|
[1224] | 257 | where $T_c$, is a parameter with dimensions of time which characterizes the force, |
---|
| 258 | $\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density, and $\rm {\bf M}$ |
---|
| 259 | represents the collected contributions of the Coriolis, hydrostatic pressure gradient, |
---|
| 260 | non-linear and viscous terms in \eqref{Eq_PE_dyn}. |
---|
[707] | 261 | |
---|
[1224] | 262 | The new force can be interpreted as a diffusion of vertically integrated volume flux divergence. |
---|
| 263 | The time evolution of $D$ is thus governed by a balance of two terms, $-g$ \textbf{A} $\eta$ |
---|
| 264 | and $g \, T_c \,$ \textbf{A} $D$, associated with a propagative regime and a diffusive regime |
---|
| 265 | in the temporal spectrum, respectively. In the diffusive regime, the EGWs no longer propagate, |
---|
| 266 | $i.e.$ they are stationary and damped. The diffusion regime applies to the modes shorter than |
---|
| 267 | $T_c$. For longer ones, the diffusion term vanishes. Hence, the temporally unresolved EGWs |
---|
| 268 | can be damped by choosing $T_c > \Delta t$. \citet{Roullet2000} demonstrate that |
---|
| 269 | (\ref{Eq_PE_flt}) can be integrated with a leap frog scheme except the additional term which |
---|
| 270 | has to be computed implicitly. This is not surprising since the use of a large time step has a |
---|
| 271 | necessarily numerical cost. Two gains arise in comparison with the previous formulations. |
---|
| 272 | Firstly, the damping of EGWs can be quantified through the magnitude of the additional term. |
---|
| 273 | Secondly, the numerical scheme does not need any tuning. Numerical stability is ensured as |
---|
| 274 | soon as $T_c > \Delta t$. |
---|
[707] | 275 | |
---|
[1224] | 276 | When the variations of free surface elevation are small compared to the thickness of the first |
---|
| 277 | model layer, the free surface equation (\ref{Eq_PE_ssh}) can be linearized. As emphasized |
---|
| 278 | by \citet{Roullet2000} the linearization of (\ref{Eq_PE_ssh}) has consequences on the |
---|
| 279 | conservation of salt in the model. With the nonlinear free surface equation, the time evolution |
---|
| 280 | of the total salt content is |
---|
[707] | 281 | \begin{equation} \label{Eq_PE_salt_content} |
---|
[1224] | 282 | \frac{\partial }{\partial t}\int\limits_{D\eta } {S\;dv} |
---|
| 283 | =\int\limits_S {S\;(-\frac{\partial \eta }{\partial t}-D+P-E)\;ds} |
---|
[707] | 284 | \end{equation} |
---|
[1224] | 285 | where $S$ is the salinity, and the total salt is integrated over the whole ocean volume |
---|
| 286 | $D_\eta$ bounded by the time-dependent free surface. The right hand side (which is an |
---|
| 287 | integral over the free surface) vanishes when the nonlinear equation (\ref{Eq_PE_ssh}) |
---|
| 288 | is satisfied, so that the salt is perfectly conserved. When the free surface equation is |
---|
| 289 | linearized, \citet{Roullet2000} show that the total salt content integrated in the fixed |
---|
| 290 | volume $D$ (bounded by the surface $z=0$) is no longer conserved: |
---|
[707] | 291 | \begin{equation} \label{Eq_PE_salt_content_linear} |
---|
[1224] | 292 | \frac{\partial }{\partial t}\int\limits_D {S\;dv} |
---|
| 293 | = - \int\limits_S {S\;\frac{\partial \eta }{\partial t}ds} |
---|
[707] | 294 | \end{equation} |
---|
| 295 | |
---|
[1224] | 296 | The right hand side of (\ref{Eq_PE_salt_content_linear}) is small in equilibrium solutions |
---|
| 297 | \citep{Roullet2000}. It can be significant when the freshwater forcing is not balanced and |
---|
| 298 | the globally averaged free surface is drifting. An increase in sea surface height \textit{$\eta $} |
---|
| 299 | results in a decrease of the salinity in the fixed volume $D$. Even in that case though, |
---|
| 300 | the total salt integrated in the variable volume $D_{\eta}$ varies much less, since |
---|
| 301 | (\ref{Eq_PE_salt_content_linear}) can be rewritten as |
---|
[707] | 302 | \begin{equation} \label{Eq_PE_salt_content_corrected} |
---|
| 303 | \frac{\partial }{\partial t}\int\limits_{D\eta } {S\;dv} |
---|
| 304 | =\frac{\partial}{\partial t} \left[ \;{\int\limits_D {S\;dv} +\int\limits_S {S\eta \;ds} } \right] |
---|
| 305 | =\int\limits_S {\eta \;\frac{\partial S}{\partial t}ds} |
---|
| 306 | \end{equation} |
---|
| 307 | |
---|
[1224] | 308 | Although the total salt content is not exactly conserved with the linearized free surface, |
---|
| 309 | its variations are driven by correlations of the time variation of surface salinity with the |
---|
| 310 | sea surface height, which is a negligible term. This situation contrasts with the case of |
---|
| 311 | the rigid lid approximation (following section) in which case freshwater forcing is |
---|
| 312 | represented by a virtual salt flux, leading to a spurious source of salt at the ocean |
---|
| 313 | surface \citep{Roullet2000}. |
---|
[707] | 314 | |
---|
| 315 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 316 | % Rigid-Lid Formulation |
---|
| 317 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 318 | \subsection{Rigid-Lid formulation} |
---|
| 319 | \label{PE_rigid_lid} |
---|
| 320 | |
---|
[1224] | 321 | With the rigid lid approximation, we assume that the ocean surface ($z=0$) is a rigid lid |
---|
| 322 | on which a pressure $p_s$ is exerted. This implies that the vertical velocity at the surface |
---|
| 323 | is equal to zero. From the continuity equation \eqref{Eq_PE_continuity} and the kinematic |
---|
| 324 | condition at the bottom \eqref{Eq_PE_w_bbc} (no flux across the bottom), it can be shown |
---|
| 325 | that the vertically integrated flow $H{\rm {\bf \bar {U}}}_h$ is non-divergent (where the |
---|
| 326 | overbar indicates a vertical average over the whole water column, i.e. from $z=-H$, |
---|
| 327 | the ocean bottom, to $z=0$ , the rigid-lid). Thus, $\rm {\bf \bar {U}}_h$ can be derived |
---|
| 328 | from a volume transport streamfunction $\psi$: |
---|
[707] | 329 | \begin{equation} \label{Eq_PE_u_psi} |
---|
| 330 | \overline{\vect{U}}_h =\frac{1}{H}\left( \vect{k} \times \nabla \psi \right) |
---|
| 331 | \end{equation} |
---|
| 332 | |
---|
[1224] | 333 | As $p_s$ does not depend on depth, its horizontal gradient is obtained by forming the |
---|
| 334 | vertical average of \eqref{Eq_PE_dyn} and using \eqref{Eq_PE_u_psi}: |
---|
[707] | 335 | |
---|
| 336 | \begin{equation} \label{Eq_PE_u_barotrope} |
---|
| 337 | \frac{1}{\rho _o }\nabla _h p_s |
---|
| 338 | =\overline{\vect{M}} -\frac{\partial \overline{\vect{U}} _h }{\partial t} |
---|
| 339 | =\overline{\vect{M}} |
---|
| 340 | -\frac{1}{H} \left[ \vect{k} \times \nabla \left( \frac{\partial \psi}{\partial t} \right) \right] |
---|
| 341 | \end{equation} |
---|
| 342 | |
---|
[1224] | 343 | Here ${\rm {\bf M}} = \left( M_u,M_v \right)$ represents the collected contributions of the |
---|
| 344 | Coriolis, hydrostatic pressure gradient, nonlinear and viscous terms in \eqref{Eq_PE_dyn}. |
---|
| 345 | The time derivative of $\psi $ is the solution of an elliptic equation which is obtained from |
---|
| 346 | the vertical component of the curl of (\ref{Eq_PE_u_barotrope}): |
---|
[707] | 347 | \begin{equation} \label{Eq_PE_psi} |
---|
| 348 | \left[ {\nabla \times \left[ {\frac{1}{H} \vect{\bf k} |
---|
| 349 | \times \nabla \left( {\frac{\partial \psi }{\partial t}} \right)} \right]} \; \right]_z |
---|
| 350 | =\left[ {\nabla \times \overline{\vect{M}} } \right]_z |
---|
| 351 | \end{equation} |
---|
| 352 | |
---|
[1224] | 353 | Using the proper boundary conditions, \eqref{Eq_PE_psi} can be solved to find $\partial_t \psi$ |
---|
| 354 | and thus using \eqref{Eq_PE_u_barotrope} the horizontal surface pressure gradient. |
---|
| 355 | It should be noted that $p_s$ can be computed by taking the divergence of |
---|
| 356 | \eqref{Eq_PE_u_barotrope} and solving the resulting elliptic equation. Thus the surface |
---|
| 357 | pressure is a diagnostic quantity that can be recovered for analysis purposes. |
---|
[707] | 358 | |
---|
[1224] | 359 | A difficulty lies in the determination of the boundary condition on $\partial_t \psi$. |
---|
| 360 | The boundary condition on velocity is that there is no flow normal to a solid wall, |
---|
| 361 | $i.e.$ the coastlines are streamlines. Therefore \eqref{Eq_PE_psi} is solved with |
---|
| 362 | the following Dirichlet boundary condition: $\partial_t \psi$ is constant along each |
---|
| 363 | coastline of the same continent or of the same island. When all the coastlines are |
---|
| 364 | connected (there are no islands), the constant value of $\partial_t \psi$ along the |
---|
| 365 | coast can be arbitrarily chosen to be zero. When islands are present in the domain, |
---|
| 366 | the value of the barotropic streamfunction will generally be different for each island |
---|
| 367 | and for the continent, and will vary with respect to time. So the boundary condition is: |
---|
| 368 | $\psi=0$ along the continent and $\psi=\mu_n$ along island $n$ ($1 \leq n \leq Q$), |
---|
| 369 | where $Q$ is the number of islands present in the domain and $\mu_n$ is a time |
---|
| 370 | dependent variable. A time evolution equation of the unknown $\mu_n$ can be found |
---|
| 371 | by evaluating the circulation of the time derivative of the vertical average (barotropic) |
---|
| 372 | velocity field along a closed contour around each island. Since the circulation of a |
---|
| 373 | gradient field along a closed contour is zero, from \eqref{Eq_PE_u_barotrope} we have: |
---|
[707] | 374 | \begin{equation} \label{Eq_PE_isl_circulation} |
---|
| 375 | \oint_n {\frac{1}{H}\left[ {{\rm {\bf k}}\times \nabla \left( |
---|
| 376 | {\frac{\partial \psi }{\partial t}} \right)} \right] \cdot {\rm {\bf d}}\ell } |
---|
| 377 | = \oint_n {\overline {\rm {\bf M}} \cdot {\rm {\bf d}}\ell } |
---|
| 378 | \qquad 1 \leq n \leq Q |
---|
| 379 | \end{equation} |
---|
| 380 | |
---|
[817] | 381 | Since (\ref{Eq_PE_psi}) is linear, its solution \textit{$\psi $} can be decomposed |
---|
| 382 | as follows: |
---|
[707] | 383 | \begin{equation} \label{Eq_PE_psi_isl} |
---|
| 384 | \psi =\psi _o +\sum\limits_{n=1}^{n=Q} {\mu _n \psi _n } |
---|
| 385 | \end{equation} |
---|
[817] | 386 | where $\psi _o$ is the solution of \eqref{Eq_PE_psi} with $\psi _o=0$ long all |
---|
| 387 | the coastlines, and where $\psi _n$ is the solution of \eqref{Eq_PE_psi} with |
---|
| 388 | the right-hand side equal to $0$, and with $\psi _n =1$ long the island $n$, |
---|
| 389 | $\psi _n =0$ along the other boundaries. The function $\psi _n$ is thus |
---|
| 390 | independent of time. Introducing \eqref{Eq_PE_psi_isl} into |
---|
| 391 | \eqref{Eq_PE_isl_circulation} yields: |
---|
[707] | 392 | \begin{multline} \label{Eq_PE_psi_isl_circulation} |
---|
[817] | 393 | \left[ {\oint_n {\frac{1}{H} \left[ {{\rm {\bf k}}\times \nabla \psi _m } \right]\cdot {\rm {\bf d}}\ell } } \right]_{1\leq m\leqslant Q \atop 1\leq n\leqslant Q } |
---|
[707] | 394 | \left( {\frac{\partial \mu _n }{\partial t}} |
---|
| 395 | \right)_{1\leqslant n\leqslant Q} \\ |
---|
| 396 | =\left( {\oint_n {\left[ {\overline {\rm |
---|
| 397 | {\bf M}} -\frac{1}{H}\left[ {{\rm {\bf k}}\times \nabla \left( |
---|
| 398 | {\frac{\partial \psi _o }{\partial t}} \right)} \right]} \right]\cdot {\rm |
---|
| 399 | {\bf d}}\ell } } \right)_{1\leqslant n\leqslant Q} |
---|
| 400 | \end{multline} |
---|
| 401 | which can be rewritten as: |
---|
| 402 | \begin{equation} \label{Eq_PE_psi_isl_matrix} |
---|
| 403 | {\rm {\bf A}}\;\left( {\frac{\partial \mu _n }{\partial t}} |
---|
| 404 | \right)_{1\leqslant n\leqslant Q} ={\rm {\bf B}} |
---|
| 405 | \end{equation} |
---|
[1224] | 406 | where \textbf{A} is a $Q \times Q$ matrix and \textbf{B} is a time dependent vector. |
---|
| 407 | As \textbf{A} is independent of time, it can be calculated and inverted once. The time |
---|
| 408 | derivative of the streamfunction when islands are present is thus given by: |
---|
[707] | 409 | \begin{equation} \label{Eq_PE_psi_isl_dt} |
---|
| 410 | \frac{\partial \psi }{\partial t}=\frac{\partial \psi _o }{\partial |
---|
| 411 | t}+\sum\limits_{n=1}^{n=Q} {{\rm {\bf A}}^{-1}{\rm {\bf B}}\;\psi _n } |
---|
| 412 | \end{equation} |
---|
| 413 | |
---|
| 414 | |
---|
| 415 | |
---|
| 416 | % ================================================================ |
---|
| 417 | % Curvilinear z-coordinate System |
---|
| 418 | % ================================================================ |
---|
| 419 | \section{Curvilinear \textit{z-}coordinate System} |
---|
| 420 | \label{PE_zco} |
---|
| 421 | |
---|
| 422 | |
---|
| 423 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 424 | % Tensorial Formalism |
---|
| 425 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 426 | \subsection{Tensorial Formalism} |
---|
| 427 | \label{PE_tensorial} |
---|
| 428 | |
---|
[1224] | 429 | In many ocean circulation problems, the flow field has regions of enhanced dynamics |
---|
| 430 | ($i.e.$ surface layers, western boundary currents, equatorial currents, or ocean fronts). |
---|
| 431 | The representation of such dynamical processes can be improved by specifically increasing |
---|
| 432 | the model resolution in these regions. As well, it may be convenient to use a lateral |
---|
| 433 | boundary-following coordinate system to better represent coastal dynamics. Moreover, |
---|
| 434 | the common geographical coordinate system has a singular point at the North Pole that |
---|
| 435 | cannot be easily treated in a global model without filtering. A solution consists of introducing |
---|
| 436 | an appropriate coordinate transformation that shifts the singular point onto land |
---|
| 437 | \citep{MadecImb1996, Murray1996}. As a consequence, it is important to solve the primitive |
---|
| 438 | equations in various curvilinear coordinate systems. An efficient way of introducing an |
---|
| 439 | appropriate coordinate transform can be found when using a tensorial formalism. |
---|
| 440 | This formalism is suited to any multidimensional curvilinear coordinate system. Ocean |
---|
| 441 | modellers mainly use three-dimensional orthogonal grids on the sphere (spherical earth |
---|
| 442 | approximation), with preservation of the local vertical. Here we give the simplified equations |
---|
| 443 | for this particular case. The general case is detailed by \citet{Eiseman1980} in their survey |
---|
| 444 | of the conservation laws of fluid dynamics. |
---|
[707] | 445 | |
---|
[1224] | 446 | Let (\textbf{i},\textbf{j},\textbf{k}) be a set of orthogonal curvilinear coordinates on the sphere |
---|
| 447 | associated with the positively oriented orthogonal set of unit vectors (\textbf{i},\textbf{j},\textbf{k}) |
---|
| 448 | linked to the earth such that \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are |
---|
| 449 | two vectors orthogonal to \textbf{k}, $i.e.$ along geopotential surfaces (Fig.\ref{Fig_referential}). |
---|
| 450 | Let $(\lambda,\varphi,z)$ be the geographical coordinate system in which a position is defined |
---|
| 451 | by the latitude $\varphi(i,j)$, the longitude $\lambda(i,j)$ and the distance from the centre of |
---|
| 452 | the earth $a+z(k)$ where $a$ is the earth's radius and $z$ the altitude above a reference sea |
---|
| 453 | level (Fig.\ref{Fig_referential}). The local deformation of the curvilinear coordinate system is |
---|
| 454 | given by $e_1$, $e_2$ and $e_3$, the three scale factors: |
---|
[707] | 455 | \begin{equation} \label{Eq_scale_factors} |
---|
| 456 | \begin{aligned} |
---|
| 457 | e_1 &=\left( {a+z} \right)\;\left[ {\left( {\frac{\partial \lambda |
---|
| 458 | }{\partial i}\cos \varphi } \right)^2+\left( {\frac{\partial \varphi |
---|
| 459 | }{\partial i}} \right)^2} \right]^{1/2} \\ |
---|
| 460 | e_2 &=\left( {a+z} \right)\;\left[ {\left( {\frac{\partial \lambda |
---|
| 461 | }{\partial j}\cos \varphi } \right)^2+\left( {\frac{\partial \varphi |
---|
| 462 | }{\partial j}} \right)^2} \right]^{1/2} \\ |
---|
| 463 | e_3 &=\left( {\frac{\partial z}{\partial k}} \right) \\ |
---|
| 464 | \end{aligned} |
---|
| 465 | \end{equation} |
---|
| 466 | |
---|
| 467 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 468 | \begin{figure}[!tb] \label{Fig_referential} \begin{center} |
---|
[998] | 469 | \includegraphics[width=0.60\textwidth]{./TexFiles/Figures/Fig_I_earth_referential.pdf} |
---|
[1224] | 470 | \caption{the geographical coordinate system $(\lambda,\varphi,z)$ and the curvilinear |
---|
| 471 | coordinate system (\textbf{i},\textbf{j},\textbf{k}). } |
---|
[707] | 472 | \end{center} \end{figure} |
---|
| 473 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 474 | |
---|
[1224] | 475 | Since the ocean depth is far smaller than the earth's radius, $a+z$, can be replaced by |
---|
| 476 | $a$ in (\ref{Eq_scale_factors}) (thin-shell approximation). The resulting horizontal scale |
---|
| 477 | factors $e_1$, $e_2$ are independent of $k$ while the vertical scale factor is a single |
---|
| 478 | function of $k$ as \textbf{k} is parallel to \textbf{z}. The scalar and vector operators that |
---|
| 479 | appear in the primitive equations (Eqs. \eqref{Eq_PE_dyn} to \eqref{Eq_PE_eos}) can |
---|
| 480 | be written in the tensorial form, invariant in any orthogonal horizontal curvilinear coordinate |
---|
| 481 | system transformation: |
---|
[707] | 482 | \begin{subequations} \label{Eq_PE_discrete_operators} |
---|
| 483 | \begin{equation} \label{Eq_PE_grad} |
---|
| 484 | \nabla q=\frac{1}{e_1 }\frac{\partial q}{\partial i}\;{\rm {\bf |
---|
| 485 | i}}+\frac{1}{e_2 }\frac{\partial q}{\partial j}\;{\rm {\bf j}}+\frac{1}{e_3 |
---|
| 486 | }\frac{\partial q}{\partial k}\;{\rm {\bf k}} \\ |
---|
| 487 | \end{equation} |
---|
| 488 | \begin{equation} \label{Eq_PE_div} |
---|
| 489 | \nabla \cdot {\rm {\bf A}} |
---|
| 490 | = \frac{1}{e_1 \; e_2} \left[ |
---|
| 491 | \frac{\partial \left(e_2 \; a_1\right)}{\partial i } |
---|
| 492 | +\frac{\partial \left(e_1 \; a_2\right)}{\partial j } \right] |
---|
| 493 | + \frac{1}{e_3} \left[ \frac{\partial a_3}{\partial k } \right] |
---|
| 494 | \end{equation} |
---|
| 495 | \begin{equation} \label{Eq_PE_curl} |
---|
| 496 | \begin{split} |
---|
| 497 | \nabla \times \vect{A} = |
---|
| 498 | \left[ {\frac{1}{e_2 }\frac{\partial a_3}{\partial j} |
---|
| 499 | -\frac{1}{e_3 }\frac{\partial a_2 }{\partial k}} \right] \; \vect{i} |
---|
| 500 | &+\left[ {\frac{1}{e_3 }\frac{\partial a_1 }{\partial k} |
---|
| 501 | -\frac{1}{e_1 }\frac{\partial a_3 }{\partial i}} \right] \; \vect{j} \\ |
---|
| 502 | &+\frac{1}{e_1 e_2 } \left[ {\frac{\partial \left( {e_2 a_2 } \right)}{\partial i} |
---|
| 503 | -\frac{\partial \left( {e_1 a_1 } \right)}{\partial j}} \right] \; \vect{k} |
---|
| 504 | \end{split} |
---|
| 505 | \end{equation} |
---|
| 506 | \begin{equation} \label{Eq_PE_lap} |
---|
| 507 | \Delta q = \nabla \cdot \left( \nabla q \right) |
---|
| 508 | \end{equation} |
---|
| 509 | \begin{equation} \label{Eq_PE_lap_vector} |
---|
| 510 | \Delta {\rm {\bf A}} = |
---|
| 511 | \nabla \left( \nabla \cdot {\rm {\bf A}} \right) |
---|
| 512 | - \nabla \times \left( \nabla \times {\rm {\bf A}} \right) |
---|
| 513 | \end{equation} |
---|
| 514 | \end{subequations} |
---|
[817] | 515 | where $q$ is a scalar quantity and ${\rm {\bf A}}=(a_1,a_2,a_3)$ a vector in the $(i,j,k)$ coordinate system. |
---|
[707] | 516 | |
---|
| 517 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 518 | % Continuous Model Equations |
---|
| 519 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 520 | \subsection{Continuous Model Equations} |
---|
| 521 | \label{PE_zco_Eq} |
---|
| 522 | |
---|
[1224] | 523 | In order to express the Primitive Equations in tensorial formalism, it is necessary to compute |
---|
| 524 | the horizontal component of the non-linear and viscous terms of the equation using |
---|
| 525 | \eqref{Eq_PE_grad}) to \eqref{Eq_PE_lap_vector}. |
---|
| 526 | Let us set $\vect U=(u,v,w)={\vect{U}}_h +w\;\vect{k}$, the velocity in the $(i,j,k)$ coordinate |
---|
| 527 | system and define the relative vorticity $\zeta$ and the divergence of the horizontal velocity |
---|
| 528 | field $\chi$, by: |
---|
[707] | 529 | \begin{equation} \label{Eq_PE_curl_Uh} |
---|
| 530 | \zeta =\frac{1}{e_1 e_2 }\left[ {\frac{\partial \left( {e_2 \,v} |
---|
| 531 | \right)}{\partial i}-\frac{\partial \left( {e_1 \,u} \right)}{\partial j}} |
---|
| 532 | \right] |
---|
| 533 | \end{equation} |
---|
| 534 | \begin{equation} \label{Eq_PE_div_Uh} |
---|
| 535 | \chi =\frac{1}{e_1 e_2 }\left[ {\frac{\partial \left( {e_2 \,u} |
---|
| 536 | \right)}{\partial i}+\frac{\partial \left( {e_1 \,v} \right)}{\partial j}} |
---|
| 537 | \right] |
---|
| 538 | \end{equation} |
---|
| 539 | |
---|
[1224] | 540 | Using the fact that the horizontal scale factors $e_1$ and $e_2$ are independent of $k$ |
---|
| 541 | and that $e_3$ is a function of the single variable $k$, the nonlinear term of |
---|
| 542 | \eqref{Eq_PE_dyn} can be transformed as follows: |
---|
[707] | 543 | \begin{flalign*} |
---|
| 544 | &\left[ {\left( { \nabla \times {\rm {\bf U}} } \right) \times {\rm {\bf U}} |
---|
| 545 | +\frac{1}{2} \nabla \left( {{\rm {\bf U}}^2} \right)} \right]_h & |
---|
| 546 | \end{flalign*} |
---|
| 547 | \begin{flalign*} |
---|
[817] | 548 | &\qquad=\left( {{\begin{array}{*{20}c} |
---|
[707] | 549 | {\left[ { \frac{1}{e_3} \frac{\partial u }{\partial k} |
---|
| 550 | -\frac{1}{e_1} \frac{\partial w }{\partial i} } \right] w - \zeta \; v } \\ |
---|
| 551 | {\zeta \; u - \left[ { \frac{1}{e_2} \frac{\partial w}{\partial j} |
---|
| 552 | -\frac{1}{e_3} \frac{\partial v}{\partial k} } \right] \ w} \\ |
---|
| 553 | \end{array} }} \right) |
---|
| 554 | +\frac{1}{2} \left( {{\begin{array}{*{20}c} |
---|
| 555 | { \frac{1}{e_1} \frac{\partial \left( u^2+v^2+w^2 \right)}{\partial i}} \hfill \\ |
---|
| 556 | { \frac{1}{e_2} \frac{\partial \left( u^2+v^2+w^2 \right)}{\partial j}} \hfill \\ |
---|
| 557 | \end{array} }} \right) & |
---|
| 558 | \end{flalign*} |
---|
| 559 | \begin{flalign*} |
---|
[817] | 560 | & \qquad =\left( {{ \begin{array}{*{20}c} |
---|
[707] | 561 | {-\zeta \; v} \hfill \\ |
---|
| 562 | { \zeta \; u} \hfill \\ |
---|
| 563 | \end{array} }} \right) |
---|
| 564 | +\frac{1}{2}\left( {{ \begin{array}{*{20}c} |
---|
| 565 | {\frac{1}{e_1 }\frac{\partial \left( {u^2+v^2} \right)}{\partial i}} \hfill \\ |
---|
| 566 | {\frac{1}{e_2 }\frac{\partial \left( {u^2+v^2} \right)}{\partial j}} \hfill \\ |
---|
| 567 | \end{array} }} \right) |
---|
| 568 | +\frac{1}{e_3 }\left( {{ \begin{array}{*{20}c} |
---|
| 569 | { w \; \frac{\partial u}{\partial k}} \\ |
---|
| 570 | { w \; \frac{\partial v}{\partial k}} \\ |
---|
| 571 | \end{array} }} \right) |
---|
| 572 | -\left( {{ \begin{array}{*{20}c} |
---|
| 573 | {\frac{w}{e_1}\frac{\partial w}{\partial i} |
---|
| 574 | -\frac{1}{2e_1}\frac{\partial w^2}{\partial i}} \hfill \\ |
---|
| 575 | {\frac{w}{e_2}\frac{\partial w}{\partial j} |
---|
| 576 | -\frac{1}{2e_2}\frac{\partial w^2}{\partial j}} \hfill \\ |
---|
| 577 | \end{array} }} \right) & |
---|
| 578 | \end{flalign*} |
---|
| 579 | |
---|
[1224] | 580 | The last term of the right hand side is obviously zero, and thus the nonlinear term of |
---|
| 581 | \eqref{Eq_PE_dyn} is written in the $(i,j,k)$ coordinate system: |
---|
[707] | 582 | \begin{equation} \label{Eq_PE_vector_form} |
---|
| 583 | \left[ {\left( { \nabla \times {\rm {\bf U}} } \right) \times {\rm {\bf U}} |
---|
| 584 | +\frac{1}{2} \nabla \left( {{\rm {\bf U}}^2} \right)} \right]_h |
---|
| 585 | =\zeta |
---|
| 586 | \;{\rm {\bf k}}\times {\rm {\bf U}}_h +\frac{1}{2}\nabla _h \left( {{\rm |
---|
| 587 | {\bf U}}_h^2 } \right)+\frac{1}{e_3 }w\frac{\partial {\rm {\bf U}}_h |
---|
| 588 | }{\partial k} |
---|
| 589 | \end{equation} |
---|
| 590 | |
---|
[1224] | 591 | This is the so-called \textit{vector invariant form} of the momentum advection term. |
---|
| 592 | For some purposes, it can be advantageous to write this term in the so-called flux form, |
---|
| 593 | $i.e.$ to write it as the divergence of fluxes. For example, the first component of |
---|
| 594 | \eqref{Eq_PE_vector_form} (the $i$-component) is transformed as follows: |
---|
[707] | 595 | \begin{flalign*} |
---|
| 596 | &{ \begin{array}{*{20}l} |
---|
| 597 | \left[ {\left( {\nabla \times \vect{U}} \right)\times \vect{U} |
---|
[817] | 598 | +\frac{1}{2}\nabla \left( {\vect{U}}^2 \right)} \right]_i % \\ |
---|
| 599 | %\\ |
---|
[707] | 600 | = - \zeta \;v |
---|
| 601 | + \frac{1}{2\;e_1 } \frac{\partial \left( {u^2+v^2} \right)}{\partial i} |
---|
| 602 | + \frac{1}{e_3}w \ \frac{\partial u}{\partial k} \\ |
---|
| 603 | \\ |
---|
[817] | 604 | \qquad =\frac{1}{e_1 \; e_2} \left( -v\frac{\partial \left( {e_2 \,v} \right)}{\partial i} |
---|
[707] | 605 | +v\frac{\partial \left( {e_1 \,u} \right)}{\partial j} \right) |
---|
| 606 | +\frac{1}{e_1 e_2 }\left( +e_2 \; u\frac{\partial u}{\partial i} |
---|
| 607 | +e_2 \; v\frac{\partial v}{\partial i} \right) |
---|
| 608 | +\frac{1}{e_3} \left( w\;\frac{\partial u}{\partial k} \right) \\ |
---|
| 609 | \end{array} } & |
---|
| 610 | \end{flalign*} |
---|
| 611 | \begin{flalign*} |
---|
| 612 | &{ \begin{array}{*{20}l} |
---|
[817] | 613 | \qquad =\frac{1}{e_1 \; e_2} \left\{ |
---|
[707] | 614 | -\left( v^2 \frac{\partial e_2 }{\partial i} |
---|
| 615 | +e_2 \,v \frac{\partial v }{\partial i} \right) |
---|
| 616 | +\left( \frac{\partial \left( {e_1 \,u\,v} \right)}{\partial j} |
---|
| 617 | -e_1 \,u \frac{\partial v }{\partial j} \right) \right. |
---|
| 618 | \\ \left. \qquad \qquad \quad |
---|
| 619 | +\left( \frac{\partial \left( {e_2 u\,u} \right)}{\partial i} |
---|
| 620 | -u \frac{\partial \left( {e_2 u} \right)}{\partial i} \right) |
---|
| 621 | +e_2 v \frac{\partial v }{\partial i} |
---|
| 622 | \right\} |
---|
| 623 | +\frac{1}{e_3} \left( |
---|
[817] | 624 | \frac{\partial \left( {w\,u} \right) }{\partial k} |
---|
[707] | 625 | -u \frac{\partial w }{\partial k} \right) \\ |
---|
| 626 | \end{array} } & |
---|
| 627 | \end{flalign*} |
---|
| 628 | \begin{flalign*} |
---|
| 629 | &{ \begin{array}{*{20}l} |
---|
[817] | 630 | \qquad =\frac{1}{e_1 \; e_2} \left( |
---|
[707] | 631 | \frac{\partial \left( {e_2 \,u\,u} \right)}{\partial i} |
---|
| 632 | + \frac{\partial \left( {e_1 \,u\,v} \right)}{\partial j} \right) |
---|
[817] | 633 | +\frac{1}{e_3 } \frac{\partial \left( {w\,u } \right)}{\partial k} |
---|
[707] | 634 | \\ \qquad \qquad \quad |
---|
| 635 | +\frac{1}{e_1 e_2 } \left( |
---|
| 636 | -u \left( \frac{\partial \left( {e_1 v } \right)}{\partial j} |
---|
| 637 | -v\,\frac{\partial e_1 }{\partial j} \right) |
---|
| 638 | -u \frac{\partial \left( {e_2 u } \right)}{\partial i} |
---|
| 639 | \right) |
---|
| 640 | -\frac{1}{e_3 } \frac{\partial w}{\partial k} u |
---|
| 641 | +\frac{1}{e_1 e_2 }\left( -v^2\frac{\partial e_2 }{\partial i} \right) |
---|
| 642 | \end{array} } & |
---|
| 643 | \end{flalign*} |
---|
| 644 | \begin{flalign*} |
---|
| 645 | &{ \begin{array}{*{20}l} |
---|
[817] | 646 | \qquad = \nabla \cdot \left( {{\rm {\bf U}}\,u} \right) |
---|
| 647 | - \left( \nabla \cdot {\rm {\bf U}} \right) \ u |
---|
[707] | 648 | +\frac{1}{e_1 e_2 }\left( |
---|
| 649 | -v^2 \frac{\partial e_2 }{\partial i} |
---|
| 650 | +uv \, \frac{\partial e_1 }{\partial j} \right) \\ |
---|
| 651 | \end{array} } & |
---|
| 652 | \end{flalign*} |
---|
| 653 | as $\nabla \cdot {\rm {\bf U}}\;=0$ (incompressibility) it comes: |
---|
| 654 | \begin{flalign*} |
---|
| 655 | &{ \begin{array}{*{20}l} |
---|
[817] | 656 | \qquad = \nabla \cdot \left( {{\rm {\bf U}}\,u} \right) |
---|
[707] | 657 | + \frac{1}{e_1 e_2 } \left( v \; \frac{\partial e_2}{\partial i} |
---|
| 658 | -u \; \frac{\partial e_1}{\partial j} \right) \left( -v \right) |
---|
| 659 | \end{array} } & |
---|
| 660 | \end{flalign*} |
---|
| 661 | |
---|
[817] | 662 | The flux form of the momentum advection term is therefore given by: |
---|
[707] | 663 | \begin{multline} \label{Eq_PE_flux_form} |
---|
| 664 | \left[ |
---|
| 665 | \left( {\nabla \times {\rm {\bf U}}} \right) \times {\rm {\bf U}} |
---|
| 666 | +\frac{1}{2} \nabla \left( {{\rm {\bf U}}^2} \right) |
---|
| 667 | \right]_h |
---|
[817] | 668 | \\ |
---|
| 669 | = \nabla \cdot \left( {{\begin{array}{*{20}c} {\rm {\bf U}} \, u \hfill \\ |
---|
[707] | 670 | {\rm {\bf U}} \, v \hfill \\ |
---|
| 671 | \end{array} }} |
---|
| 672 | \right) |
---|
| 673 | +\frac{1}{e_1 e_2 } \left( |
---|
| 674 | v\frac{\partial e_2}{\partial i} |
---|
| 675 | -u\frac{\partial e_1}{\partial j} |
---|
| 676 | \right) {\rm {\bf k}} \times {\rm {\bf U}}_h |
---|
| 677 | \end{multline} |
---|
| 678 | |
---|
[1224] | 679 | The flux form has two terms, the first one is expressed as the divergence of momentum |
---|
| 680 | fluxes (hence the flux form name given to this formulation) and the second one is due to |
---|
| 681 | the curvilinear nature of the coordinate system used. The latter is called the \emph{metric} |
---|
| 682 | term and can be viewed as a modification of the Coriolis parameter: |
---|
[707] | 683 | \begin{equation} \label{Eq_PE_cor+metric} |
---|
| 684 | f \to f + \frac{1}{e_1 \; e_2} \left( v \frac{\partial e_2}{\partial i} |
---|
| 685 | -u \frac{\partial e_1}{\partial j} \right) |
---|
| 686 | \end{equation} |
---|
| 687 | |
---|
[1224] | 688 | Note that in the case of geographical coordinate, $i.e.$ when $(i,j) \to (\lambda ,\varphi )$ |
---|
| 689 | and $(e_1 ,e_2) \to (a \,\cos \varphi ,a)$, we recover the commonly used modification of |
---|
| 690 | the Coriolis parameter $f \to f+(u/a) \tan \varphi$. |
---|
[707] | 691 | |
---|
[817] | 692 | To sum up, the equations solved by the ocean model can be written in the following tensorial formalism: |
---|
[707] | 693 | |
---|
[817] | 694 | \vspace{+10pt} |
---|
| 695 | $\bullet$ \textit{momentum equations} : |
---|
| 696 | |
---|
| 697 | vector invariant form : |
---|
[707] | 698 | \begin{subequations} \label{Eq_PE_dyn_vect} |
---|
[817] | 699 | \begin{multline} \label{Eq_PE_dyn_vect_u} |
---|
[707] | 700 | \frac{\partial u}{\partial t}= |
---|
[817] | 701 | + \left( {\zeta +f} \right)\,v |
---|
| 702 | - \frac{1}{2\,e_1} \frac{\partial}{\partial i} \left( u^2+v^2 \right) |
---|
| 703 | - \frac{1}{e_3} w \frac{\partial u}{\partial k} \\ |
---|
| 704 | - \frac{1}{e_1} \frac{\partial}{\partial i} \left( \frac{p_s+p_h }{\rho _o} \right) |
---|
| 705 | + D_u^{\vect{U}} + F_u^{\vect{U}} |
---|
| 706 | \end{multline} |
---|
| 707 | \begin{multline} \label{Eq_PE_dyn_vect_v} |
---|
[707] | 708 | \frac{\partial v}{\partial t}= |
---|
[817] | 709 | - \left( {\zeta +f} \right)\,u |
---|
| 710 | - \frac{1}{2\,e_2 }\frac{\partial }{\partial j}\left( u^2+v^2 \right) - \frac{1}{e_3 }w\frac{\partial v}{\partial k} \\ |
---|
| 711 | - \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho _o} \right) |
---|
| 712 | + D_v^{\vect{U}} + F_v^{\vect{U}} |
---|
| 713 | \end{multline} |
---|
[707] | 714 | \end{subequations} |
---|
| 715 | |
---|
[817] | 716 | flux form: |
---|
| 717 | \begin{subequations} \label{Eq_PE_dyn_flux} |
---|
| 718 | \begin{multline} \label{Eq_PE_dyn_flux_u} |
---|
| 719 | \frac{\partial u}{\partial t}= |
---|
| 720 | + \left( { f + \frac{1}{e_1 \; e_2} |
---|
| 721 | \left( v \frac{\partial e_2}{\partial i} |
---|
| 722 | -u \frac{\partial e_1}{\partial j} \right)} \right) \, v \\ |
---|
| 723 | - \frac{1}{e_1 \; e_2} \left( |
---|
| 724 | \frac{\partial \left( {e_2 \,u\,u} \right)}{\partial i} |
---|
| 725 | + \frac{\partial \left( {e_1 \,v\,u} \right)}{\partial j} \right) |
---|
[994] | 726 | - \frac{1}{e_3 }\frac{\partial \left( { w\,u} \right)}{\partial k} \\ |
---|
[817] | 727 | - \frac{1}{e_1 }\frac{\partial}{\partial i}\left( \frac{p_s+p_h }{\rho _o} \right) |
---|
| 728 | + D_u^{\vect{U}} + F_u^{\vect{U}} |
---|
| 729 | \end{multline} |
---|
| 730 | \begin{multline} \label{Eq_PE_dyn_flux_v} |
---|
| 731 | \frac{\partial v}{\partial t}= |
---|
| 732 | - \left( { f + \frac{1}{e_1 \; e_2} |
---|
| 733 | \left( v \frac{\partial e_2}{\partial i} |
---|
| 734 | -u \frac{\partial e_1}{\partial j} \right)} \right) \, u \\ |
---|
| 735 | \frac{1}{e_1 \; e_2} \left( |
---|
| 736 | \frac{\partial \left( {e_2 \,u\,v} \right)}{\partial i} |
---|
[994] | 737 | + \frac{\partial \left( {e_1 \,v\,v} \right)}{\partial j} \right) |
---|
| 738 | - \frac{1}{e_3 } \frac{\partial \left( { w\,v} \right)}{\partial k} \\ |
---|
[817] | 739 | - \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho _o} \right) |
---|
| 740 | + D_v^{\vect{U}} + F_v^{\vect{U}} |
---|
| 741 | \end{multline} |
---|
| 742 | \end{subequations} |
---|
[1224] | 743 | where $\zeta$ is given by \eqref{Eq_PE_curl_Uh} and the surface pressure gradient formulation |
---|
| 744 | depends on the one of the free surface: |
---|
[707] | 745 | |
---|
| 746 | $*$ free surface formulation |
---|
| 747 | \begin{equation}\label{Eq_PE_dyn_sco} |
---|
| 748 | \frac{1}{\rho _o }\nabla _h p_s =\left( {{\begin{array}{*{20}c} |
---|
| 749 | {\frac{g}{\;e_1 }\frac{\partial \eta }{\partial i}} \hfill \\ |
---|
| 750 | {\frac{g}{\;e_2 }\frac{\partial \eta }{\partial j}} \hfill \\ |
---|
| 751 | \end{array} }} \right) |
---|
| 752 | \qquad \text{where $\eta$ is solution of \eqref{Eq_PE_ssh} } |
---|
| 753 | \end{equation} |
---|
| 754 | |
---|
| 755 | $*$ rigid-lid approximation |
---|
| 756 | \begin{equation}\label{Eq_PE_dyn_zco} |
---|
| 757 | \frac{1}{\rho _o }\nabla _h p_s =\left( {{\begin{array}{*{20}c} |
---|
| 758 | {\overline M _u +\frac{1}{H\;e_2 }\frac{\partial }{\partial j}\left( |
---|
| 759 | {\frac{\partial \psi }{\partial t}} \right)} \\ |
---|
| 760 | {\overline M _v -\frac{1}{H\;e_1 }\frac{\partial }{\partial i}\left( |
---|
| 761 | {\frac{\partial \psi }{\partial t}} \right)} \\ |
---|
| 762 | \end{array} }} \right) |
---|
| 763 | \end{equation} |
---|
| 764 | where ${\vect{M}}= \left( M_u,M_v \right)$ represents the collected contributions of nonlinear, |
---|
[1224] | 765 | viscosity and hydrostatic pressure gradient terms in \eqref{Eq_PE_dyn_vect} and the overbar |
---|
| 766 | indicates a vertical average over the whole water column ($i.e.$ from $z=-H$, the ocean bottom, |
---|
| 767 | to $z=0$, the rigid-lid), and where the time derivative of $\psi$ is the solution of an elliptic equation: |
---|
[707] | 768 | \begin{multline} \label{Eq_psi_total} |
---|
| 769 | \frac{\partial }{\partial i}\left[ {\frac{e_2 }{H\,e_1}\frac{\partial}{\partial i} |
---|
| 770 | \left( {\frac{\partial \psi }{\partial t}} \right)} \right] |
---|
| 771 | +\frac{\partial }{\partial j}\left[ {\frac{e_1 }{H\,e_2}\frac{\partial }{\partial j} |
---|
| 772 | \left( {\frac{\partial \psi }{\partial t}} \right)} \right] |
---|
| 773 | = \\ |
---|
| 774 | \frac{\partial }{\partial i}\left( {e_2 \overline M _v } \right) |
---|
| 775 | - \frac{\partial }{\partial j}\left( {e_1 \overline M _u } \right) |
---|
| 776 | \end{multline} |
---|
| 777 | |
---|
| 778 | The vertical velocity and the hydrostatic pressure are diagnosed from the following equations: |
---|
| 779 | \begin{equation} \label{Eq_w_diag} |
---|
| 780 | \frac{\partial w}{\partial k}=-\chi \;e_3 |
---|
| 781 | \end{equation} |
---|
| 782 | \begin{equation} \label{Eq_hp_diag} |
---|
| 783 | \frac{\partial p_h }{\partial k}=-\rho \;g\;e_3 |
---|
| 784 | \end{equation} |
---|
| 785 | |
---|
[817] | 786 | where the divergence of the horizontal velocity, $\chi$ is given by \eqref{Eq_PE_div_Uh}. |
---|
[707] | 787 | |
---|
[817] | 788 | \vspace{+10pt} |
---|
| 789 | $\bullet$ \textit{tracer equations} : |
---|
[707] | 790 | \begin{equation} \label{Eq_S} |
---|
[817] | 791 | \frac{\partial T}{\partial t} = |
---|
| 792 | -\frac{1}{e_1 e_2 }\left[ { \frac{\partial \left( {e_2 T\,u} \right)}{\partial i} |
---|
| 793 | +\frac{\partial \left( {e_1 T\,v} \right)}{\partial j}} \right] |
---|
| 794 | -\frac{1}{e_3 }\frac{\partial \left( {T\,w} \right)}{\partial k} + D^T + F^T |
---|
[707] | 795 | \end{equation} |
---|
| 796 | \begin{equation} \label{Eq_T} |
---|
[817] | 797 | \frac{\partial S}{\partial t} = |
---|
| 798 | -\frac{1}{e_1 e_2 }\left[ {\frac{\partial \left( {e_2 S\,u} \right)}{\partial i} |
---|
| 799 | +\frac{\partial \left( {e_1 S\,v} \right)}{\partial j}} \right] |
---|
| 800 | -\frac{1}{e_3 }\frac{\partial \left( {S\,w} \right)}{\partial k} + D^S + F^S |
---|
[707] | 801 | \end{equation} |
---|
| 802 | \begin{equation} \label{Eq_rho} |
---|
| 803 | \rho =\rho \left( {T,S,z(k)} \right) |
---|
| 804 | \end{equation} |
---|
| 805 | |
---|
[1224] | 806 | The expression of \textbf{D}$^{U}$, $D^{S}$ and $D^{T}$ depends on the subgrid scale |
---|
| 807 | parameterisation used. It will be defined in \S\ref{PE_zdf}. The nature and formulation of |
---|
| 808 | ${\rm {\bf F}}^{\rm {\bf U}}$, $F^T$ and $F^S$, the surface forcing terms, are discussed |
---|
| 809 | in Chapter~\ref{SBC}. |
---|
[707] | 810 | |
---|
| 811 | \newpage |
---|
| 812 | % ================================================================ |
---|
[817] | 813 | % Curvilinear z*-coordinate System |
---|
| 814 | % ================================================================ |
---|
| 815 | \section{Curvilinear \textit{z*}-coordinate System} |
---|
| 816 | |
---|
| 817 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 818 | \begin{figure}[!b] \label{Fig_z_zstar} \begin{center} |
---|
[998] | 819 | \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_z_zstar.pdf} |
---|
[1224] | 820 | \caption{(a) $z$-coordinate in linear free-surface case ; (b) $z-$coordinate in non-linear |
---|
| 821 | free surface case (c) re-scaled height coordinate (become popular as the \textit{z*-}coordinate |
---|
| 822 | \citep{Adcroft_Campin_OM04} ).} |
---|
[817] | 823 | \end{center} \end{figure} |
---|
| 824 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 825 | |
---|
| 826 | |
---|
[1224] | 827 | In that case, the free surface equation is nonlinear, and the variations of volume are fully |
---|
| 828 | taken into account. These coordinates systems is presented in a report \citep{Levier2007} |
---|
| 829 | available on the \NEMO web site. |
---|
[817] | 830 | |
---|
| 831 | \gmcomment{ |
---|
[1224] | 832 | The \textit{z*} coordinate approach is an unapproximated, non-linear free surface implementation |
---|
| 833 | which allows one to deal with large amplitude free-surface |
---|
[817] | 834 | variations relative to the vertical resolution \citep{Adcroft_Campin_OM04}. In |
---|
| 835 | the \textit{z*} formulation, the variation of the column thickness due to sea-surface |
---|
| 836 | undulations is not concentrated in the surface level, as in the $z$-coordinate formulation, |
---|
| 837 | but is equally distributed over the full water column. Thus vertical |
---|
| 838 | levels naturally follow sea-surface variations, with a linear attenuation with |
---|
| 839 | depth, as illustrated by figure fig.1c . Note that with a flat bottom, such as in |
---|
| 840 | fig.1c, the bottom-following $z$ coordinate and \textit{z*} are equivalent. |
---|
| 841 | The definition and modified oceanic equations for the rescaled vertical coordinate |
---|
| 842 | \textit{z*}, including the treatment of fresh-water flux at the surface, are |
---|
| 843 | detailed in Adcroft and Campin (2004). The major points are summarized |
---|
| 844 | here. The position ( \textit{z*}) and vertical discretization (\textit{z*}) are expressed as: |
---|
| 845 | |
---|
| 846 | $H + \textit{z*} = (H + z) / r$ and $\delta \textit{z*} = \delta z / r$ with $r = \frac{H+\eta} {H}$ |
---|
| 847 | |
---|
| 848 | Since the vertical displacement of the free surface is incorporated in the vertical |
---|
[1224] | 849 | coordinate \textit{z*}, the upper and lower boundaries are at fixed \textit{z*} position, |
---|
| 850 | $\textit{z*} = 0$ and $\textit{z*} = ?H$ respectively. Also the divergence of the flow field |
---|
| 851 | is no longer zero as shown by the continuity equation: |
---|
[817] | 852 | |
---|
| 853 | $\frac{\partial r}{\partial t} = \nabla_{\textit{z*}} \cdot \left( r \; \rm{\bf U}_h \right) |
---|
| 854 | \left( r \; w\textit{*} \right) = 0 $ |
---|
| 855 | |
---|
| 856 | } |
---|
| 857 | |
---|
| 858 | |
---|
| 859 | \newpage |
---|
| 860 | % ================================================================ |
---|
[707] | 861 | % Curvilinear s-coordinate System |
---|
| 862 | % ================================================================ |
---|
| 863 | \section{Curvilinear \textit{s}-coordinate System} |
---|
| 864 | \label{PE_sco} |
---|
| 865 | |
---|
| 866 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 867 | % Introduction |
---|
| 868 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 869 | \subsection{Introduction} |
---|
| 870 | |
---|
[1224] | 871 | Several important aspects of the ocean circulation are influenced by bottom topography. |
---|
| 872 | Of course, the most important is that bottom topography determines deep ocean sub-basins, |
---|
| 873 | barriers, sills and channels that strongly constrain the path of water masses, but more subtle |
---|
| 874 | effects exist. For example, the topographic $\beta$-effect is usually larger than the planetary |
---|
| 875 | one along continental slopes. Topographic Rossby waves can be excited and can interact |
---|
| 876 | with the mean current. In the $z-$coordinate system presented in the previous section |
---|
| 877 | (\S\ref{PE_zco}), $z-$surfaces are geopotential surfaces. The bottom topography is |
---|
| 878 | discretised by steps. This often leads to a misrepresentation of a gradually sloping bottom |
---|
| 879 | and to large localized depth gradients associated with large localized vertical velocities. |
---|
| 880 | The response to such a velocity field often leads to numerical dispersion effects. |
---|
| 881 | One solution to strongly reduce this error is to use a partial step representation of bottom |
---|
| 882 | topography instead of a full step one \cite{Pacanowski_Gnanadesikan_MWR98}. |
---|
| 883 | Another solution is to introduce a terrain-following coordinate system (hereafter $s-$coordinate) |
---|
[707] | 884 | |
---|
[1224] | 885 | The $s$-coordinate avoids the discretisation error in the depth field since the layers of |
---|
| 886 | computation are gradually adjusted with depth to the ocean bottom. Relatively small |
---|
| 887 | topographic features as well as gentle, large-scale slopes of the sea floor in the deep |
---|
| 888 | ocean, which would be ignored in typical $z$-model applications with the largest grid |
---|
| 889 | spacing at greatest depths, can easily be represented (with relatively low vertical resolution). |
---|
| 890 | A terrain-following model (hereafter $s-$model) also facilitates the modelling of the |
---|
| 891 | boundary layer flows over a large depth range, which in the framework of the $z$-model |
---|
| 892 | would require high vertical resolution over the whole depth range. Moreover, with a |
---|
| 893 | $s$-coordinate it is possible, at least in principle, to have the bottom and the sea surface |
---|
| 894 | as the only boundaries of the domain (nomore lateral boundary condition to specify). |
---|
| 895 | Nevertheless, a $s$-coordinate also has its drawbacks. Perfectly adapted to a |
---|
| 896 | homogeneous ocean, it has strong limitations as soon as stratification is introduced. |
---|
| 897 | The main two problems come from the truncation error in the horizontal pressure |
---|
| 898 | gradient and a possibly increased diapycnal diffusion. The horizontal pressure force |
---|
| 899 | in $s$-coordinate consists of two terms (see Appendix~\ref{Apdx_A}), |
---|
[707] | 900 | |
---|
| 901 | \begin{equation} \label{Eq_PE_p_sco} |
---|
| 902 | \left. {\nabla p} \right|_z =\left. {\nabla p} \right|_s -\frac{\partial |
---|
| 903 | p}{\partial s}\left. {\nabla z} \right|_s |
---|
| 904 | \end{equation} |
---|
| 905 | |
---|
[1224] | 906 | The second term in \eqref{Eq_PE_p_sco} depends on the tilt of the coordinate surface |
---|
| 907 | and introduces a truncation error that is not present in a $z$-model. In the special case |
---|
| 908 | of a $\sigma-$coordinate (i.e. a depth-normalised coordinate system $\sigma = z/H$), |
---|
| 909 | \citet{Haney1991} and \citet{Beckmann1993} have given estimates of the magnitude |
---|
| 910 | of this truncation error. It depends on topographic slope, stratification, horizontal and |
---|
| 911 | vertical resolution, the equation of state, and the finite difference scheme. This error |
---|
| 912 | limits the possible topographic slopes that a model can handle at a given horizontal |
---|
| 913 | and vertical resolution. This is a severe restriction for large-scale applications using |
---|
| 914 | realistic bottom topography. The large-scale slopes require high horizontal resolution, |
---|
| 915 | and the computational cost becomes prohibitive. This problem can be at least partially |
---|
| 916 | overcome by mixing $s$-coordinate and step-like representation of bottom topography \citep{Gerdes1993a,Gerdes1993b,Madec1996}. However, the definition of the model |
---|
| 917 | domain vertical coordinate becomes then a non-trivial thing for a realistic bottom |
---|
| 918 | topography: a envelope topography is defined in $s$-coordinate on which a full or |
---|
| 919 | partial step bottom topography is then applied in order to adjust the model depth to |
---|
| 920 | the observed one (see \S\ref{DOM_zgr}. |
---|
[707] | 921 | |
---|
[1224] | 922 | For numerical reasons a minimum of diffusion is required along the coordinate surfaces |
---|
| 923 | of any finite difference model. It causes spurious diapycnal mixing when coordinate |
---|
| 924 | surfaces do not coincide with isoneutral surfaces. This is the case for a $z$-model as |
---|
| 925 | well as for a $s$-model. However, density varies more strongly on $s-$surfaces than |
---|
| 926 | on horizontal surfaces in regions of large topographic slopes, implying larger diapycnal |
---|
| 927 | diffusion in a $s$-model than in a $z$-model. Whereas such a diapycnal diffusion in a |
---|
| 928 | $z$-model tends to weaken horizontal density (pressure) gradients and thus the horizontal |
---|
| 929 | circulation, it usually reinforces these gradients in a $s$-model, creating spurious circulation. |
---|
| 930 | For example, imagine an isolated bump of topography in an ocean at rest with a horizontally |
---|
| 931 | uniform stratification. Spurious diffusion along $s$-surfaces will induce a bump of isoneutral |
---|
| 932 | surfaces over the topography, and thus will generate there a baroclinic eddy. In contrast, |
---|
| 933 | the ocean will stay at rest in a $z$-model. As for the truncation error, the problem can be reduced by introducing the terrain-following coordinate below the strongly stratified portion of the water column |
---|
| 934 | ($i.e.$ the main thermocline) \citep{Madec1996}. An alternate solution consists of rotating |
---|
| 935 | the lateral diffusive tensor to geopotential or to isoneutral surfaces (see \S\ref{PE_ldf}. |
---|
| 936 | Unfortunately, the slope of isoneutral surfaces relative to the $s$-surfaces can very large, |
---|
| 937 | strongly exceeding the stability limit of such a operator when it is discretized (see Chapter~\ref{LDF}). |
---|
[707] | 938 | |
---|
[1224] | 939 | The $s-$coordinates introduced here \citep{Lott1990,Madec1996} differ mainly in two |
---|
| 940 | aspects from similar models: it allows a representation of bottom topography with mixed |
---|
| 941 | full or partial step-like/terrain following topography ; It also offers a completely general |
---|
| 942 | transformation, $s=s(i,j,z)$ for the vertical coordinate. |
---|
[707] | 943 | |
---|
| 944 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 945 | % The s-coordinate Formulation |
---|
| 946 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 947 | \subsection{The \textit{s-}coordinate Formulation} |
---|
| 948 | |
---|
[1224] | 949 | Starting from the set of equations established in \S\ref{PE_zco} for the special case $k=z$ |
---|
| 950 | and thus $e_3=1$, we introduce an arbitrary vertical coordinate $s=s(i,j,k,t)$, which includes |
---|
| 951 | $z$-, \textit{z*}- and $\sigma-$coordinates as special cases ($s=z$, $s=\textit{z*}$, and |
---|
| 952 | $s=\sigma=z/H$ or $=z/\left(H+\eta \right)$, resp.). A formal derivation of the transformed |
---|
| 953 | equations is given in Appendix~\ref{Apdx_A}. Let us define the vertical scale factor by |
---|
| 954 | $e_3=\partial_s z$ ($e_3$ is now a function of $(i,j,k,t)$ ), and the slopes in the |
---|
| 955 | (\textbf{i},\textbf{j}) directions between $s-$ and $z-$surfaces by : |
---|
[707] | 956 | \begin{equation} \label{Eq_PE_sco_slope} |
---|
| 957 | \sigma _1 =\frac{1}{e_1 }\;\left. {\frac{\partial z}{\partial i}} \right|_s |
---|
| 958 | \quad \text{, and } \quad |
---|
| 959 | \sigma _2 =\frac{1}{e_2 }\;\left. {\frac{\partial z}{\partial j}} \right|_s |
---|
| 960 | \end{equation} |
---|
[1224] | 961 | We also introduce $\omega $, a dia-surface velocity component, defined as the velocity |
---|
| 962 | relative to the moving $s$-surfaces and normal to them: |
---|
[707] | 963 | \begin{equation} \label{Eq_PE_sco_w} |
---|
[817] | 964 | \omega = w - e_3 \, \frac{\partial s}{\partial t} - \sigma _1 \,u - \sigma _2 \,v \\ |
---|
[707] | 965 | \end{equation} |
---|
| 966 | |
---|
[817] | 967 | The equations solved by the ocean model \eqref{Eq_PE} in $s-$coordinate can be written as follows: |
---|
| 968 | |
---|
| 969 | \vspace{0.5cm} |
---|
[707] | 970 | * momentum equation: |
---|
| 971 | \begin{multline} \label{Eq_PE_sco_u} |
---|
[817] | 972 | \frac{1}{e_3} \frac{\partial \left( e_3\,u \right) }{\partial t}= |
---|
| 973 | + \left( {\zeta +f} \right)\,v |
---|
| 974 | - \frac{1}{2\,e_1} \frac{\partial}{\partial i} \left( u^2+v^2 \right) |
---|
| 975 | - \frac{1}{e_3} \omega \frac{\partial u}{\partial k} \\ |
---|
| 976 | - \frac{1}{e_1} \frac{\partial}{\partial i} \left( \frac{p_s + p_h}{\rho _o} \right) |
---|
| 977 | + g\frac{\rho }{\rho _o}\sigma _1 |
---|
| 978 | + D_u^{\vect{U}} + F_u^{\vect{U}} \quad |
---|
[707] | 979 | \end{multline} |
---|
| 980 | \begin{multline} \label{Eq_PE_sco_v} |
---|
[817] | 981 | \frac{1}{e_3} \frac{\partial \left( e_3\,v \right) }{\partial t}= |
---|
| 982 | - \left( {\zeta +f} \right)\,u |
---|
| 983 | - \frac{1}{2\,e_2 }\frac{\partial }{\partial j}\left( u^2+v^2 \right) |
---|
| 984 | - \frac{1}{e_3 } \omega \frac{\partial v}{\partial k} \\ |
---|
| 985 | - \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho _o} \right) |
---|
| 986 | + g\frac{\rho }{\rho _o }\sigma _2 |
---|
| 987 | + D_v^{\vect{U}} + F_v^{\vect{U}} \quad |
---|
[707] | 988 | \end{multline} |
---|
[1224] | 989 | where the relative vorticity, \textit{$\zeta $}, the surface pressure gradient, and the hydrostatic |
---|
| 990 | pressure have the same expressions as in $z$-coordinates although they do not represent |
---|
| 991 | exactly the same quantities. $\omega$ is provided by the continuity equation |
---|
| 992 | (see Appendix~\ref{Apdx_A}): |
---|
[707] | 993 | |
---|
[817] | 994 | \begin{equation} \label{Eq_PE_sco_continuity} |
---|
| 995 | \frac{\partial e_3}{\partial t} + e_3 \; \chi + \frac{\partial \omega }{\partial s} = 0 |
---|
| 996 | \qquad \text{with }\;\; |
---|
[707] | 997 | \chi =\frac{1}{e_1 e_2 e_3 }\left[ {\frac{\partial \left( {e_2 e_3 \,u} |
---|
| 998 | \right)}{\partial i}+\frac{\partial \left( {e_1 e_3 \,v} \right)}{\partial |
---|
| 999 | j}} \right] |
---|
| 1000 | \end{equation} |
---|
| 1001 | |
---|
[817] | 1002 | \vspace{0.5cm} |
---|
[707] | 1003 | * tracer equations: |
---|
[817] | 1004 | \begin{multline} \label{Eq_PE_sco_t} |
---|
| 1005 | \frac{1}{e_3} \frac{\partial \left( e_3\,T \right) }{\partial t}= |
---|
| 1006 | -\frac{1}{e_1 e_2 e_3 }\left[ {\frac{\partial \left( {e_2 e_3\,u\,T} \right)}{\partial i} |
---|
| 1007 | +\frac{\partial \left( {e_1 e_3\,v\,T} \right)}{\partial j}} \right] \\ |
---|
| 1008 | -\frac{1}{e_3 }\frac{\partial \left( {T\,\omega } \right)}{\partial k} + D^T + F^S \qquad |
---|
| 1009 | \end{multline} |
---|
[707] | 1010 | |
---|
[817] | 1011 | \begin{multline} \label{Eq_PE_sco_s} |
---|
| 1012 | \frac{1}{e_3} \frac{\partial \left( e_3\,S \right) }{\partial t}= |
---|
| 1013 | -\frac{1}{e_1 e_2 e_3 }\left[ {\frac{\partial \left( {e_2 e_3\,u\,S} \right)}{\partial i} |
---|
| 1014 | +\frac{\partial \left( {e_1 e_3\,v\,S} \right)}{\partial j}} \right] \\ |
---|
| 1015 | -\frac{1}{e_3 }\frac{\partial \left( {S\,\omega } \right)}{\partial k} + D^S + F^S \qquad |
---|
| 1016 | \end{multline} |
---|
[707] | 1017 | |
---|
[1224] | 1018 | The equation of state has the same expression as in $z$-coordinate, and similar expressions |
---|
| 1019 | are used for mixing and forcing terms. |
---|
[707] | 1020 | |
---|
[817] | 1021 | \gmcomment{ |
---|
[707] | 1022 | \colorbox{yellow}{ to be updated $= = >$} |
---|
| 1023 | Add a few works on z and zps and s and underlies the differences between all of them |
---|
[817] | 1024 | \colorbox{yellow}{ $< = =$ end update} } |
---|
[707] | 1025 | |
---|
| 1026 | \newpage |
---|
| 1027 | % ================================================================ |
---|
| 1028 | % Subgrid Scale Physics |
---|
| 1029 | % ================================================================ |
---|
| 1030 | \section{Subgrid Scale Physics} |
---|
| 1031 | \label{PE_zdf_ldf} |
---|
| 1032 | |
---|
| 1033 | The primitive equations describe the behaviour of a geophysical fluid at |
---|
| 1034 | space and time scales larger than a few kilometres in the horizontal, a few |
---|
| 1035 | meters in the vertical and a few minutes. They are usually solved at larger |
---|
[817] | 1036 | scales: the specified grid spacing and time step of the numerical model. The |
---|
[707] | 1037 | effects of smaller scale motions (coming from the advective terms in the |
---|
| 1038 | Navier-Stokes equations) must be represented entirely in terms of |
---|
| 1039 | large-scale patterns to close the equations. These effects appear in the |
---|
[817] | 1040 | equations as the divergence of turbulent fluxes ($i.e.$ fluxes associated with |
---|
[707] | 1041 | the mean correlation of small scale perturbations). Assuming a turbulent |
---|
| 1042 | closure hypothesis is equivalent to choose a formulation for these fluxes. |
---|
| 1043 | It is usually called the subgrid scale physics. It must be emphasized that |
---|
| 1044 | this is the weakest part of the primitive equations, but also one of the |
---|
[1224] | 1045 | most important for long-term simulations as small scale processes \textit{in fine} |
---|
| 1046 | balance the surface input of kinetic energy and heat. |
---|
[707] | 1047 | |
---|
| 1048 | The control exerted by gravity on the flow induces a strong anisotropy |
---|
[1224] | 1049 | between the lateral and vertical motions. Therefore subgrid-scale physics |
---|
| 1050 | \textbf{D}$^{\vect{U}}$, $D^{S}$ and $D^{T}$ in \eqref{Eq_PE_dyn}, |
---|
| 1051 | \eqref{Eq_PE_tra_T} and \eqref{Eq_PE_tra_S} are divided into a lateral part |
---|
| 1052 | \textbf{D}$^{l \vect{U}}$, $D^{lS}$ and $D^{lT}$ and a vertical part |
---|
| 1053 | \textbf{D}$^{vU}$, $D^{vS}$ and $D^{vT}$. The formulation of these terms |
---|
| 1054 | and their underlying physics are briefly discussed in the next two subsections. |
---|
[707] | 1055 | |
---|
| 1056 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 1057 | % Vertical Subgrid Scale Physics |
---|
| 1058 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 1059 | \subsection{Vertical Subgrid Scale Physics} |
---|
| 1060 | \label{PE_zdf} |
---|
| 1061 | |
---|
| 1062 | The model resolution is always larger than the scale at which the major |
---|
[817] | 1063 | sources of vertical turbulence occur (shear instability, internal wave |
---|
[707] | 1064 | breaking...). Turbulent motions are thus never explicitly solved, even |
---|
| 1065 | partially, but always parameterized. The vertical turbulent fluxes are |
---|
| 1066 | assumed to depend linearly on the gradients of large-scale quantities (for |
---|
[1224] | 1067 | example, the turbulent heat flux is given by $\overline{T'w'}=-A^{vT} \partial_z \overline T$, |
---|
| 1068 | where $A^{vT}$ is an eddy coefficient). This formulation is |
---|
[707] | 1069 | analogous to that of molecular diffusion and dissipation. This is quite |
---|
| 1070 | clearly a necessary compromise: considering only the molecular viscosity |
---|
| 1071 | acting on large scale severely underestimates the role of turbulent |
---|
| 1072 | diffusion and dissipation, while an accurate consideration of the details of |
---|
| 1073 | turbulent motions is simply impractical. The resulting vertical momentum and |
---|
| 1074 | tracer diffusive operators are of second order: |
---|
| 1075 | \begin{equation} \label{Eq_PE_zdf} |
---|
| 1076 | \begin{split} |
---|
[1224] | 1077 | {\vect{D}}^{v \vect{U}} &=\frac{\partial }{\partial z}\left( {A^{vm}\frac{\partial {\vect{U}}_h }{\partial z}} \right) \ , \\ |
---|
| 1078 | D^{vT} &= \frac{\partial }{\partial z}\left( {A^{vT}\frac{\partial T}{\partial z}} \right) \ , |
---|
[707] | 1079 | \quad |
---|
| 1080 | D^{vS}=\frac{\partial }{\partial z}\left( {A^{vT}\frac{\partial S}{\partial z}} \right) |
---|
| 1081 | \end{split} |
---|
| 1082 | \end{equation} |
---|
[1224] | 1083 | where $A^{vm}$ and $A^{vT}$ are the vertical eddy viscosity and diffusivity coefficients, |
---|
| 1084 | respectively. At the sea surface and at the bottom, turbulent fluxes of momentum, heat |
---|
| 1085 | and salt must be specified (see Chap.~\ref{SBC} and \ref{ZDF} and \S\ref{TRA_bbl}). |
---|
| 1086 | All the vertical physics is embedded in the specification of the eddy coefficients. |
---|
| 1087 | They can be assumed to be either constant, or function of the local fluid properties |
---|
| 1088 | ($e.g.$ Richardson number, Brunt-Vais\"{a}l\"{a} frequency...), or computed from a |
---|
| 1089 | turbulent closure model. The choices available in \NEMO are discussed in \S\ref{ZDF}). |
---|
[707] | 1090 | |
---|
| 1091 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 1092 | % Lateral Diffusive and Viscous Operators Formulation |
---|
| 1093 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 1094 | \subsection{Lateral Diffusive and Viscous Operators Formulation} |
---|
| 1095 | \label{PE_ldf} |
---|
| 1096 | |
---|
| 1097 | Lateral turbulence can be roughly divided into a mesoscale turbulence |
---|
[817] | 1098 | associated with eddies (which can be solved explicitly if the resolution is |
---|
| 1099 | sufficient since their underlying physics are included in the primitive |
---|
| 1100 | equations), and a sub mesoscale turbulence which is never explicitly solved |
---|
[707] | 1101 | even partially, but always parameterized. The formulation of lateral eddy |
---|
| 1102 | fluxes depends on whether the mesoscale is below or above the grid-spacing |
---|
[817] | 1103 | ($i.e.$ the model is eddy-resolving or not). |
---|
[707] | 1104 | |
---|
[817] | 1105 | In non-eddy-resolving configurations, the closure is similar to that used |
---|
[707] | 1106 | for the vertical physics. The lateral turbulent fluxes are assumed to depend |
---|
| 1107 | linearly on the lateral gradients of large-scale quantities. The resulting |
---|
| 1108 | lateral diffusive and dissipative operators are of second order. |
---|
| 1109 | Observations show that lateral mixing induced by mesoscale turbulence tends |
---|
[1224] | 1110 | to be along isopycnal surfaces (or more precisely neutral surfaces \cite{McDougall1987}) |
---|
| 1111 | rather than across them. |
---|
[817] | 1112 | As the slope of neutral surfaces is small in the ocean, a common |
---|
[707] | 1113 | approximation is to assume that the `lateral' direction is the horizontal, |
---|
[817] | 1114 | $i.e.$ the lateral mixing is performed along geopotential surfaces. This leads |
---|
[707] | 1115 | to a geopotential second order operator for lateral subgrid scale physics. |
---|
| 1116 | This assumption can be relaxed: the eddy-induced turbulent fluxes can be |
---|
| 1117 | better approached by assuming that they depend linearly on the gradients of |
---|
[817] | 1118 | large-scale quantities computed along neutral surfaces. In such a case, |
---|
[707] | 1119 | the diffusive operator is an isoneutral second order operator and it has |
---|
| 1120 | components in the three space directions. However, both horizontal and |
---|
[817] | 1121 | isoneutral operators have no effect on mean ($i.e.$ large scale) potential |
---|
[707] | 1122 | energy whereas potential energy is a main source of turbulence (through |
---|
| 1123 | baroclinic instabilities). \citet{Gent1990} have proposed a |
---|
[1224] | 1124 | parameterisation of mesoscale eddy-induced turbulence which associates an |
---|
[707] | 1125 | eddy-induced velocity to the isoneutral diffusion. Its mean effect is to |
---|
| 1126 | reduce the mean potential energy of the ocean. This leads to a formulation |
---|
| 1127 | of lateral subgrid-scale physics made up of an isoneutral second order |
---|
| 1128 | operator and an eddy induced advective part. In all these lateral diffusive |
---|
| 1129 | formulations, the specification of the lateral eddy coefficients remains the |
---|
[817] | 1130 | problematic point as there is no really satisfactory formulation of these |
---|
[707] | 1131 | coefficients as a function of large-scale features. |
---|
| 1132 | |
---|
| 1133 | In eddy-resolving configurations, a second order operator can be used, but |
---|
| 1134 | usually a more scale selective one (biharmonic operator) is preferred as the |
---|
| 1135 | grid-spacing is usually not small enough compared to the scale of the |
---|
| 1136 | eddies. The role devoted to the subgrid-scale physics is to dissipate the |
---|
| 1137 | energy that cascades toward the grid scale and thus ensures the stability of |
---|
[817] | 1138 | the model while not interfering with the solved mesoscale activity. Another approach |
---|
| 1139 | is becoming more and more popular: instead of specifying explicitly a sub-grid scale |
---|
[1224] | 1140 | term in the momentum and tracer time evolution equations, one uses a advective |
---|
| 1141 | scheme which is diffusive enough to maintain the model stability. It must be emphasised |
---|
[817] | 1142 | that then, all the sub-grid scale physics is in this case include in the formulation of the |
---|
| 1143 | advection scheme. |
---|
[707] | 1144 | |
---|
[1224] | 1145 | All these parameterisations of subgrid scale physics present advantages and |
---|
[817] | 1146 | drawbacks. There are not all available in \NEMO. In the $z$-coordinate |
---|
| 1147 | formulation, five options are offered for active tracers (temperature and |
---|
[707] | 1148 | salinity): second order geopotential operator, second order isoneutral |
---|
[1224] | 1149 | operator, \citet{Gent1990} parameterisation, fourth order |
---|
| 1150 | geopotential operator, and various slightly diffusive advection schemes. |
---|
| 1151 | The same options are available for momentum, except |
---|
| 1152 | \citet{Gent1990} parameterisation which only involves tracers. In the |
---|
[817] | 1153 | $s$-coordinate formulation, additional options are offered for tracers: second |
---|
[707] | 1154 | order operator acting along $s-$surfaces, and for momentum: fourth order |
---|
| 1155 | operator acting along $s-$surfaces (see \S\ref{LDF}). |
---|
| 1156 | |
---|
| 1157 | \subsubsection{lateral second order tracer diffusive operator} |
---|
| 1158 | |
---|
[817] | 1159 | The lateral second order tracer diffusive operator is defined by (see Appendix~\ref{Apdx_B}): |
---|
[707] | 1160 | \begin{equation} \label{Eq_PE_iso_tensor} |
---|
| 1161 | D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad |
---|
| 1162 | \mbox{with}\quad \;\;\Re =\left( {{\begin{array}{*{20}c} |
---|
| 1163 | 1 \hfill & 0 \hfill & {-r_1 } \hfill \\ |
---|
| 1164 | 0 \hfill & 1 \hfill & {-r_2 } \hfill \\ |
---|
| 1165 | {-r_1 } \hfill & {-r_2 } \hfill & {r_1 ^2+r_2 ^2} \hfill \\ |
---|
| 1166 | \end{array} }} \right) |
---|
| 1167 | \end{equation} |
---|
| 1168 | where $r_1 \;\mbox{and}\;r_2 $ are the slopes between the surface along |
---|
[817] | 1169 | which the diffusive operator acts and the model level ($e. g.$ $z$- or |
---|
| 1170 | $s$-surfaces). Note that the formulation \eqref{Eq_PE_iso_tensor} is exact for the |
---|
| 1171 | rotation between geopotential and $s$-surfaces, while it is only an approximation |
---|
| 1172 | for the rotation between isoneutral and $z$- or $s$-surfaces. Indeed, in the latter |
---|
[1224] | 1173 | case, two assumptions are made to simplify \eqref{Eq_PE_iso_tensor} \citep{Cox1987}. |
---|
| 1174 | First, the horizontal contribution of the dianeutral mixing is neglected since the ratio |
---|
| 1175 | between iso and dia-neutral diffusive coefficients is known to be several orders of |
---|
| 1176 | magnitude smaller than unity. Second, the two isoneutral directions of diffusion are |
---|
| 1177 | assumed to be independent since the slopes are generally less than $10^{-2}$ in the |
---|
| 1178 | ocean (see Appendix~\ref{Apdx_B}). |
---|
[707] | 1179 | |
---|
| 1180 | For \textit{geopotential} diffusion, $r_1$ and $r_2 $ are the slopes between the |
---|
[1224] | 1181 | geopotential and computational surfaces: in $z$-coordinates they are zero |
---|
| 1182 | ($r_1 = r_2 = 0$) while in $s$-coordinate (including $\textit{z*}$ case) they are |
---|
| 1183 | equal to $\sigma _1$ and $\sigma _2$, respectively (see \eqref{Eq_PE_sco_slope} ). |
---|
[707] | 1184 | |
---|
[1224] | 1185 | For \textit{isoneutral} diffusion $r_1$ and $r_2$ are the slopes between the isoneutral |
---|
| 1186 | and computational surfaces. Therefore, they have a same expression in $z$- and $s$-coordinates: |
---|
[707] | 1187 | \begin{equation} \label{Eq_PE_iso_slopes} |
---|
| 1188 | r_1 =\frac{e_3 }{e_1 } \left( {\frac{\partial \rho }{\partial i}} \right) |
---|
| 1189 | \left( {\frac{\partial \rho }{\partial k}} \right)^{-1} \ , \quad |
---|
| 1190 | r_1 =\frac{e_3 }{e_1 } \left( {\frac{\partial \rho }{\partial i}} \right) |
---|
| 1191 | \left( {\frac{\partial \rho }{\partial k}} \right)^{-1} |
---|
| 1192 | \end{equation} |
---|
| 1193 | |
---|
[1224] | 1194 | When the \textit{Eddy Induced Velocity} parametrisation (eiv) \citep{Gent1990} is used, |
---|
| 1195 | an additional tracer advection is introduced in combination with the isoneutral diffusion of tracers: |
---|
[707] | 1196 | \begin{equation} \label{Eq_PE_iso+eiv} |
---|
| 1197 | D^{lT}=\nabla \cdot \left( {A^{lT}\;\Re \;\nabla T} \right) |
---|
| 1198 | +\nabla \cdot \left( {{\vect{U}}^\ast \,T} \right) |
---|
| 1199 | \end{equation} |
---|
[1224] | 1200 | where ${\vect{U}}^\ast =\left( {u^\ast ,v^\ast ,w^\ast } \right)$ is a non-divergent, |
---|
| 1201 | eddy-induced transport velocity. This velocity field is defined by: |
---|
[707] | 1202 | \begin{equation} \label{Eq_PE_eiv} |
---|
| 1203 | \begin{split} |
---|
| 1204 | u^\ast &= +\frac{1}{e_3 }\frac{\partial }{\partial k}\left[ {A^{eiv}\;\tilde{r}_1 } \right] \\ |
---|
| 1205 | v^\ast &= +\frac{1}{e_3 }\frac{\partial }{\partial k}\left[ {A^{eiv}\;\tilde{r}_2 } \right] \\ |
---|
| 1206 | w^\ast &= -\frac{1}{e_1 e_2 }\left[ |
---|
| 1207 | \frac{\partial }{\partial i}\left( {A^{eiv}\;e_2\,\tilde{r}_1 } \right) |
---|
| 1208 | +\frac{\partial }{\partial j}\left( {A^{eiv}\;e_1\,\tilde{r}_2 } \right) \right] |
---|
| 1209 | \end{split} |
---|
| 1210 | \end{equation} |
---|
[1224] | 1211 | where $A^{eiv}$ is the eddy induced velocity coefficient (or equivalently the isoneutral |
---|
| 1212 | thickness diffusivity coefficient), and $\tilde{r}_1$ and $\tilde{r}_2$ are the slopes |
---|
| 1213 | between isoneutral and \emph{geopotential} surfaces and thus depends on the coordinate |
---|
| 1214 | considered: |
---|
[707] | 1215 | \begin{align} \label{Eq_PE_slopes_eiv} |
---|
| 1216 | \tilde{r}_n = \begin{cases} |
---|
[817] | 1217 | r_n & \text{in $z$-coordinate} \\ |
---|
| 1218 | r_n + \sigma_n & \text{in \textit{z*} and $s$-coordinates} |
---|
[707] | 1219 | \end{cases} |
---|
| 1220 | \quad \text{where } n=1,2 |
---|
| 1221 | \end{align} |
---|
| 1222 | |
---|
[1224] | 1223 | The normal component of the eddy induced velocity is zero at all the boundaries. |
---|
| 1224 | This can be achieved in a model by tapering either the eddy coefficient or the slopes |
---|
| 1225 | to zero in the vicinity of the boundaries. The latter strategy is used in \NEMO (cf. Chap.~\ref{LDF}). |
---|
[707] | 1226 | |
---|
| 1227 | \subsubsection{lateral fourth order tracer diffusive operator} |
---|
| 1228 | |
---|
| 1229 | The lateral fourth order tracer diffusive operator is defined by: |
---|
| 1230 | \begin{equation} \label{Eq_PE_bilapT} |
---|
| 1231 | D^{lT}=\Delta \left( {A^{lT}\;\Delta T} \right) |
---|
| 1232 | \qquad \text{where} \ D^{lT}=\Delta \left( {A^{lT}\;\Delta T} \right) |
---|
| 1233 | \end{equation} |
---|
| 1234 | |
---|
[1224] | 1235 | It is the second order operator given by \eqref{Eq_PE_iso_tensor} applied twice with |
---|
| 1236 | the eddy diffusion coefficient correctly placed. |
---|
[707] | 1237 | |
---|
[817] | 1238 | |
---|
[707] | 1239 | \subsubsection{lateral second order momentum diffusive operator} |
---|
| 1240 | |
---|
[817] | 1241 | The second order momentum diffusive operator along $z$- or $s$-surfaces is found by |
---|
| 1242 | applying \eqref{Eq_PE_lap_vector} to the horizontal velocity vector (see Appendix~\ref{Apdx_B}): |
---|
[707] | 1243 | \begin{equation} \label{Eq_PE_lapU} |
---|
| 1244 | \begin{split} |
---|
| 1245 | {\rm {\bf D}}^{l{\rm {\bf U}}} |
---|
| 1246 | &= \quad \ \nabla _h \left( {A^{lm}\chi } \right) |
---|
| 1247 | \ - \ \nabla _h \times \left( {A^{lm}\,\zeta \;{\rm {\bf k}}} \right) \\ |
---|
| 1248 | &= \left( \begin{aligned} |
---|
| 1249 | \frac{1}{e_1 } \frac{\partial \left( A^{lm} \chi \right)}{\partial i} |
---|
| 1250 | &-\frac{1}{e_2 e_3}\frac{\partial \left( {A^{lm} \;e_3 \zeta} \right)}{\partial j} \\ |
---|
| 1251 | \frac{1}{e_2 }\frac{\partial \left( {A^{lm} \chi } \right)}{\partial j} |
---|
| 1252 | &+\frac{1}{e_1 e_3}\frac{\partial \left( {A^{lm} \;e_3 \zeta} \right)}{\partial i} |
---|
| 1253 | \end{aligned} \right) |
---|
| 1254 | \end{split} |
---|
| 1255 | \end{equation} |
---|
| 1256 | |
---|
| 1257 | Such a formulation ensures a complete separation between the vorticity and |
---|
[817] | 1258 | horizontal divergence fields (see Appendix~\ref{Apdx_C}). Unfortunately, it is not |
---|
[707] | 1259 | available for geopotential diffusion in $s-$coordinates and for isoneutral |
---|
[1224] | 1260 | diffusion in both $z$- and $s$-coordinates ($i.e.$ when a rotation is required). |
---|
| 1261 | In these two cases, the $u$ and $v-$fields are considered as independent scalar |
---|
| 1262 | fields, so that the diffusive operator is given by: |
---|
[707] | 1263 | \begin{equation} \label{Eq_PE_lapU_iso} |
---|
| 1264 | \begin{split} |
---|
| 1265 | D_u^{l{\rm {\bf U}}} &= \nabla .\left( {\Re \;\nabla u} \right) \\ |
---|
| 1266 | D_v^{l{\rm {\bf U}}} &= \nabla .\left( {\Re \;\nabla v} \right) |
---|
| 1267 | \end{split} |
---|
| 1268 | \end{equation} |
---|
[1224] | 1269 | where $\Re$ is given by \eqref{Eq_PE_iso_tensor}. It is the same expression as |
---|
| 1270 | those used for diffusive operator on tracers. It must be emphasised that such a |
---|
| 1271 | formulation is only exact in a Cartesian coordinate system, $i.e.$ on a $f-$ or |
---|
| 1272 | $\beta-$plane, not on the sphere. It is also a very good approximation in vicinity |
---|
| 1273 | of the Equator in a geographical coordinate system \citep{Lengaigne_al_JGR03}. |
---|
[707] | 1274 | |
---|
| 1275 | \subsubsection{lateral fourth order momentum diffusive operator} |
---|
| 1276 | |
---|
[1224] | 1277 | As for tracers, the fourth order momentum diffusive operator along $z$ or $s$-surfaces |
---|
| 1278 | is a re-entering second order operator \eqref{Eq_PE_lapU} or \eqref{Eq_PE_lapU} |
---|
| 1279 | with the eddy viscosity coefficient correctly placed: |
---|
[707] | 1280 | |
---|
[817] | 1281 | geopotential diffusion in $z$-coordinate: |
---|
[707] | 1282 | \begin{equation} \label{Eq_PE_bilapU} |
---|
| 1283 | \begin{split} |
---|
| 1284 | {\rm {\bf D}}^{l{\rm {\bf U}}} &=\nabla _h \left\{ {\;\nabla _h {\rm {\bf |
---|
| 1285 | .}}\left[ {A^{lm}\,\nabla _h \left( \chi \right)} \right]\;} |
---|
| 1286 | \right\}\; \\ |
---|
| 1287 | &+\nabla _h \times \left\{ {\;{\rm {\bf k}}\cdot \nabla \times |
---|
| 1288 | \left[ {A^{lm}\,\nabla _h \times \left( {\zeta \;{\rm {\bf k}}} \right)} |
---|
| 1289 | \right]\;} \right\} |
---|
| 1290 | \end{split} |
---|
| 1291 | \end{equation} |
---|
| 1292 | |
---|
[817] | 1293 | \gmcomment{ change the position of the coefficient, both here and in the code} |
---|
| 1294 | |
---|
| 1295 | geopotential diffusion in $s$-coordinate: |
---|
[707] | 1296 | \begin{equation} \label{Eq_bilapU_iso} |
---|
| 1297 | \left\{ \begin{aligned} |
---|
| 1298 | D_u^{l{\rm {\bf U}}} =\Delta \left( {A^{lm}\;\Delta u} \right) \\ |
---|
| 1299 | D_v^{l{\rm {\bf U}}} =\Delta \left( {A^{lm}\;\Delta v} \right) |
---|
| 1300 | \end{aligned} \right. |
---|
| 1301 | \quad \text{where} \quad |
---|
| 1302 | \Delta \left( \bullet \right) = \nabla \cdot \left( \Re \nabla(\bullet) \right) |
---|
| 1303 | \end{equation} |
---|
| 1304 | |
---|