- Timestamp:
- 2018-11-21T17:59:55+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics.tex
r9414 r10354 26 26 27 27 28 The ocean is a fluid that can be described to a good approximation by the primitive 29 equations, $i.e.$ the Navier-Stokes equations along with a nonlinear equation of 30 state which couples the two active tracers (temperature and salinity) to the fluid 31 velocity, plus the following additional assumptions made from scale considerations: 32 33 \textit{(1) spherical earth approximation: }the geopotential surfaces are assumed to 34 be spheres so that gravity (local vertical) is parallel to the earth's radius 35 36 \textit{(2) thin-shell approximation: }the ocean depth is neglected compared to the earth's radius 37 38 \textit{(3) turbulent closure hypothesis: }the turbulent fluxes (which represent the effect 39 of small scale processes on the large-scale) are expressed in terms of large-scale features 40 41 \textit{(4) Boussinesq hypothesis:} density variations are neglected except in their 42 contribution to the buoyancy force 43 44 \textit{(5) Hydrostatic hypothesis: }the vertical momentum equation is reduced to a 45 balance between the vertical pressure gradient and the buoyancy force (this removes 46 convective processes from the initial Navier-Stokes equations and so convective processes 47 must be parameterized instead) 48 49 \textit{(6) Incompressibility hypothesis: }the three dimensional divergence of the velocity 50 vector is assumed to be zero. 51 52 Because the gravitational force is so dominant in the equations of large-scale motions, 53 it is useful to choose an orthogonal set of unit vectors (\textbf{i},\textbf{j},\textbf{k}) linked 54 to the earth such that \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are two 55 vectors orthogonal to \textbf{k}, $i.e.$ tangent to the geopotential surfaces. Let us define 56 the following variables: \textbf{U} the vector velocity, $\textbf{U}=\textbf{U}_h + w\, \textbf{k}$ 28 The ocean is a fluid that can be described to a good approximation by the primitive equations, 29 $i.e.$ the Navier-Stokes equations along with a nonlinear equation of state which 30 couples the two active tracers (temperature and salinity) to the fluid velocity, 31 plus the following additional assumptions made from scale considerations: 32 33 \textit{(1) spherical earth approximation:} the geopotential surfaces are assumed to be spheres so that 34 gravity (local vertical) is parallel to the earth's radius 35 36 \textit{(2) thin-shell approximation:} the ocean depth is neglected compared to the earth's radius 37 38 \textit{(3) turbulent closure hypothesis:} the turbulent fluxes 39 (which represent the effect of small scale processes on the large-scale) are expressed in terms of 40 large-scale features 41 42 \textit{(4) Boussinesq hypothesis:} density variations are neglected except in their contribution to 43 the buoyancy force 44 45 \textit{(5) Hydrostatic hypothesis:} the vertical momentum equation is reduced to a balance between 46 the vertical pressure gradient and the buoyancy force 47 (this removes convective processes from the initial Navier-Stokes equations and so 48 convective processes must be parameterized instead) 49 50 \textit{(6) Incompressibility hypothesis:} the three dimensional divergence of the velocity vector is assumed to 51 be zero. 52 53 Because the gravitational force is so dominant in the equations of large-scale motions, 54 it is useful to choose an orthogonal set of unit vectors (\textbf{i},\textbf{j},\textbf{k}) linked to 55 the earth such that \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are two vectors orthogonal to 56 \textbf{k}, $i.e.$ tangent to the geopotential surfaces. 57 Let us define the following variables: \textbf{U} the vector velocity, $\textbf{U}=\textbf{U}_h + w\, \textbf{k}$ 57 58 (the subscript $h$ denotes the local horizontal vector, $i.e.$ over the (\textbf{i},\textbf{j}) plane), 58 $T$ the potential temperature, $S$ the salinity, \textit{$\rho $} the \textit{in situ} density. 59 The vector invariant form of the primitive equations in the (\textbf{i},\textbf{j},\textbf{k}) 60 vector system provides the following six equations (namely the momentum balance, the 61 hydrostatic equilibrium, the incompressibility equation, the heat and salt conservation 62 equations and an equation of state):59 $T$ the potential temperature, $S$ the salinity, \textit{$\rho $} the \textit{in situ} density. 60 The vector invariant form of the primitive equations in the (\textbf{i},\textbf{j},\textbf{k}) vector system 61 provides the following six equations 62 (namely the momentum balance, the hydrostatic equilibrium, the incompressibility equation, 63 the heat and salt conservation equations and an equation of state): 63 64 \begin{subequations} \label{eq:PE} 64 65 \begin{equation} \label{eq:PE_dyn} … … 85 86 \end{equation} 86 87 \end{subequations} 87 where $\nabla$ is the generalised derivative vector operator in $(\bf i,\bf j, \bf k)$ directions, 88 $t$ is the time, $z$ is the vertical coordinate, $\rho $ is the \textit{in situ} density given by 89 the equation of state (\autoref{eq:PE_eos}), $\rho_o$ is a reference density, $p$ the pressure, 90 $f=2 \bf \Omega \cdot \bf k$ is the Coriolis acceleration (where $\bf \Omega$ is the Earth's 91 angular velocity vector), and $g$ is the gravitational acceleration. 92 ${\rm {\bf D}}^{\rm {\bf U}}$, $D^T$ and $D^S$ are the parameterisations of small-scale 93 physics for momentum, temperature and salinity, and ${\rm {\bf F}}^{\rm {\bf U}}$, $F^T$ 94 and $F^S$ surface forcing terms. Their nature and formulation are discussed in 95 \autoref{sec:PE_zdf_ldf} and \autoref{subsec:PE_boundary_condition}. 96 97 . 88 where $\nabla$ is the generalised derivative vector operator in $(\bf i,\bf j, \bf k)$ directions, $t$ is the time, 89 $z$ is the vertical coordinate, $\rho $ is the \textit{in situ} density given by the equation of state 90 (\autoref{eq:PE_eos}), $\rho_o$ is a reference density, $p$ the pressure, 91 $f=2 \bf \Omega \cdot \bf k$ is the Coriolis acceleration 92 (where $\bf \Omega$ is the Earth's angular velocity vector), and $g$ is the gravitational acceleration. 93 ${\rm {\bf D}}^{\rm {\bf U}}$, $D^T$ and $D^S$ are the parameterisations of small-scale physics for momentum, 94 temperature and salinity, and ${\rm {\bf F}}^{\rm {\bf U}}$, $F^T$ and $F^S$ surface forcing terms. 95 Their nature and formulation are discussed in \autoref{sec:PE_zdf_ldf} and \autoref{subsec:PE_boundary_condition}. 96 97 98 98 99 99 % ------------------------------------------------------------------------------------------------------------- … … 103 103 \label{subsec:PE_boundary_condition} 104 104 105 An ocean is bounded by complex coastlines, bottom topography at its base and an air-sea 106 or ice-sea interface at its top. These boundaries can be defined by two surfaces, $z=-H(i,j)$ 107 and $z=\eta(i,j,k,t)$, where $H$ is the depth of the ocean bottom and $\eta$ is the height 108 of the sea surface. Both $H$ and $\eta$ are usually referenced to a given surface, $z=0$, 109 chosen as a mean sea surface (\autoref{fig:ocean_bc}). Through these two boundaries, 110 the ocean can exchange fluxes of heat, fresh water, salt, and momentum with the solid earth, 111 the continental margins, the sea ice and the atmosphere. However, some of these fluxes are 112 so weak that even on climatic time scales of thousands of years they can be neglected. 113 In the following, we briefly review the fluxes exchanged at the interfaces between the ocean 114 and the other components of the earth system. 105 An ocean is bounded by complex coastlines, bottom topography at its base and 106 an air-sea or ice-sea interface at its top. 107 These boundaries can be defined by two surfaces, $z=-H(i,j)$ and $z=\eta(i,j,k,t)$, 108 where $H$ is the depth of the ocean bottom and $\eta$ is the height of the sea surface. 109 Both $H$ and $\eta$ are usually referenced to a given surface, $z=0$, chosen as a mean sea surface 110 (\autoref{fig:ocean_bc}). 111 Through these two boundaries, the ocean can exchange fluxes of heat, fresh water, salt, and momentum with 112 the solid earth, the continental margins, the sea ice and the atmosphere. 113 However, some of these fluxes are so weak that even on climatic time scales of thousands of years 114 they can be neglected. 115 In the following, we briefly review the fluxes exchanged at the interfaces between the ocean and 116 the other components of the earth system. 115 117 116 118 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 117 \begin{figure}[!ht] \begin{center} 118 \includegraphics[width=0.90\textwidth]{Fig_I_ocean_bc} 119 \caption{ \protect\label{fig:ocean_bc} 120 The ocean is bounded by two surfaces, $z=-H(i,j)$ and $z=\eta(i,j,t)$, where $H$ 121 is the depth of the sea floor and $\eta$ the height of the sea surface. 122 Both $H$ and $\eta$ are referenced to $z=0$.} 123 \end{center} \end{figure} 119 \begin{figure}[!ht] 120 \begin{center} 121 \includegraphics[width=0.90\textwidth]{Fig_I_ocean_bc} 122 \caption{ \protect\label{fig:ocean_bc} 123 The ocean is bounded by two surfaces, $z=-H(i,j)$ and $z=\eta(i,j,t)$, 124 where $H$ is the depth of the sea floor and $\eta$ the height of the sea surface. 125 Both $H$ and $\eta$ are referenced to $z=0$. 126 } 127 \end{center} 128 \end{figure} 124 129 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 125 130 126 131 127 132 \begin{description} 128 \item[Land - ocean interface:] the major flux between continental margins and the ocean is 129 a mass exchange of fresh water through river runoff. Such an exchange modifies the sea 130 surface salinity especially in the vicinity of major river mouths. It can be neglected for short 131 range integrations but has to be taken into account for long term integrations as it influences 132 the characteristics of water masses formed (especially at high latitudes). It is required in order 133 to close the water cycle of the climate system. It is usually specified as a fresh water flux at 134 the air-sea interface in the vicinity of river mouths. 135 \item[Solid earth - ocean interface:] heat and salt fluxes through the sea floor are small, 136 except in special areas of little extent. They are usually neglected in the model 137 \footnote{In fact, it has been shown that the heat flux associated with the solid Earth cooling 138 ($i.e.$the geothermal heating) is not negligible for the thermohaline circulation of the world 139 ocean (see \autoref{subsec:TRA_bbc}).}. 140 The boundary condition is thus set to no flux of heat and salt across solid boundaries. 141 For momentum, the situation is different. There is no flow across solid boundaries, 142 $i.e.$ the velocity normal to the ocean bottom and coastlines is zero (in other words, 143 the bottom velocity is parallel to solid boundaries). This kinematic boundary condition 144 can be expressed as: 145 \begin{equation} \label{eq:PE_w_bbc} 146 w = -{\rm {\bf U}}_h \cdot \nabla _h \left( H \right) 147 \end{equation} 148 In addition, the ocean exchanges momentum with the earth through frictional processes. 149 Such momentum transfer occurs at small scales in a boundary layer. It must be parameterized 150 in terms of turbulent fluxes using bottom and/or lateral boundary conditions. Its specification 151 depends on the nature of the physical parameterisation used for ${\rm {\bf D}}^{\rm {\bf U}}$ 152 in \autoref{eq:PE_dyn}. It is discussed in \autoref{eq:PE_zdf}.% and Chap. III.6 to 9. 153 \item[Atmosphere - ocean interface:] the kinematic surface condition plus the mass flux 154 of fresh water PE (the precipitation minus evaporation budget) leads to: 155 \begin{equation} \label{eq:PE_w_sbc} 156 w = \frac{\partial \eta }{\partial t} 157 + \left. {{\rm {\bf U}}_h } \right|_{z=\eta } \cdot \nabla _h \left( \eta \right) 133 \item[Land - ocean interface:] 134 the major flux between continental margins and the ocean is a mass exchange of fresh water through river runoff. 135 Such an exchange modifies the sea surface salinity especially in the vicinity of major river mouths. 136 It can be neglected for short range integrations but has to be taken into account for long term integrations as 137 it influences the characteristics of water masses formed (especially at high latitudes). 138 It is required in order to close the water cycle of the climate system. 139 It is usually specified as a fresh water flux at the air-sea interface in the vicinity of river mouths. 140 \item[Solid earth - ocean interface:] 141 heat and salt fluxes through the sea floor are small, except in special areas of little extent. 142 They are usually neglected in the model \footnote{ 143 In fact, it has been shown that the heat flux associated with the solid Earth cooling 144 ($i.e.$the geothermal heating) is not negligible for the thermohaline circulation of the world ocean 145 (see \autoref{subsec:TRA_bbc}). 146 }. 147 The boundary condition is thus set to no flux of heat and salt across solid boundaries. 148 For momentum, the situation is different. There is no flow across solid boundaries, 149 $i.e.$ the velocity normal to the ocean bottom and coastlines is zero (in other words, 150 the bottom velocity is parallel to solid boundaries). This kinematic boundary condition 151 can be expressed as: 152 \begin{equation} \label{eq:PE_w_bbc} 153 w = -{\rm {\bf U}}_h \cdot \nabla _h \left( H \right) 154 \end{equation} 155 In addition, the ocean exchanges momentum with the earth through frictional processes. 156 Such momentum transfer occurs at small scales in a boundary layer. 157 It must be parameterized in terms of turbulent fluxes using bottom and/or lateral boundary conditions. 158 Its specification depends on the nature of the physical parameterisation used for 159 ${\rm {\bf D}}^{\rm {\bf U}}$ in \autoref{eq:PE_dyn}. 160 It is discussed in \autoref{eq:PE_zdf}.% and Chap. III.6 to 9. 161 \item[Atmosphere - ocean interface:] 162 the kinematic surface condition plus the mass flux of fresh water PE (the precipitation minus evaporation budget) 163 leads to: 164 \begin{equation} \label{eq:PE_w_sbc} 165 w = \frac{\partial \eta }{\partial t} 166 + \left. {{\rm {\bf U}}_h } \right|_{z=\eta } \cdot \nabla _h \left( \eta \right) 158 167 + P-E 159 \end{equation}160 The dynamic boundary condition, neglecting the surface tension (which removes capillary 161 waves from the system) leads to the continuity of pressure across the interface $z=\eta$. 162 The atmosphere and ocean also exchange horizontal momentum (wind stress), and heat.163 \item[Sea ice - ocean interface:] the ocean and sea ice exchange heat, salt, fresh water164 and momentum. The sea surface temperature is constrained to be at the freezing point 165 at the interface. Sea ice salinity is very low ($\sim4-6 \,psu$) compared to those of the 166 ocean ($\sim34 \,psu$). The cycle of freezing/melting is associated with fresh water and 167 salt fluxes that cannot be neglected.168 \end{equation} 169 The dynamic boundary condition, neglecting the surface tension (which removes capillary waves from the system) 170 leads to the continuity of pressure across the interface $z=\eta$. 171 The atmosphere and ocean also exchange horizontal momentum (wind stress), and heat. 172 \item[Sea ice - ocean interface:] 173 the ocean and sea ice exchange heat, salt, fresh water and momentum. 174 The sea surface temperature is constrained to be at the freezing point at the interface. 175 Sea ice salinity is very low ($\sim4-6 \,psu$) compared to those of the ocean ($\sim34 \,psu$). 176 The cycle of freezing/melting is associated with fresh water and salt fluxes that cannot be neglected. 168 177 \end{description} 169 178 … … 184 193 \label{subsec:PE_p_formulation} 185 194 186 The total pressure at a given depth $z$ is composed of a surface pressure $p_s$ at a 187 reference geopotential surface ($z=0$) and a hydrostatic pressure $p_h$ such that: 188 $p(i,j,k,t)=p_s(i,j,t)+p_h(i,j,k,t)$. The latter is computed by integrating (\autoref{eq:PE_hydrostatic}), 189 assuming that pressure in decibars can be approximated by depth in meters in (\autoref{eq:PE_eos}). 195 The total pressure at a given depth $z$ is composed of a surface pressure $p_s$ at 196 a reference geopotential surface ($z=0$) and a hydrostatic pressure $p_h$ such that: 197 $p(i,j,k,t)=p_s(i,j,t)+p_h(i,j,k,t)$. 198 The latter is computed by integrating (\autoref{eq:PE_hydrostatic}), 199 assuming that pressure in decibars can be approximated by depth in meters in (\autoref{eq:PE_eos}). 190 200 The hydrostatic pressure is then given by: 191 201 \begin{equation} \label{eq:PE_pressure} … … 193 203 = \int_{\varsigma =z}^{\varsigma =0} {g\;\rho \left( {T,S,\varsigma} \right)\;d\varsigma } 194 204 \end{equation} 195 Two strategies can be considered for the surface pressure term: $(a)$ introduce of a 196 new variable $\eta$, the free-surface elevation, for which a prognostic equation can be 197 established and solved; $(b)$ assume that the ocean surface is a rigid lid, on which the 198 pressure (or its horizontal gradient) can be diagnosed. When the former strategy is used, 199 one solution of the free-surface elevation consists of the excitation of external gravity waves. 200 The flow is barotropic and the surface moves up and down with gravity as the restoring force. 201 The phase speed of such waves is high (some hundreds of metres per second) so that 202 the time step would have to be very short if they were present in the model. The latter 203 strategy filters out these waves since the rigid lid approximation implies $\eta=0$, $i.e.$ 204 the sea surface is the surface $z=0$. This well known approximation increases the surface 205 wave speed to infinity and modifies certain other longwave dynamics ($e.g.$ barotropic 206 Rossby or planetary waves). The rigid-lid hypothesis is an obsolescent feature in modern 207 OGCMs. It has been available until the release 3.1 of \NEMO, and it has been removed 208 in release 3.2 and followings. Only the free surface formulation is now described in the 209 this document (see the next sub-section). 205 Two strategies can be considered for the surface pressure term: 206 $(a)$ introduce of a new variable $\eta$, the free-surface elevation, 207 for which a prognostic equation can be established and solved; 208 $(b)$ assume that the ocean surface is a rigid lid, 209 on which the pressure (or its horizontal gradient) can be diagnosed. 210 When the former strategy is used, one solution of the free-surface elevation consists of 211 the excitation of external gravity waves. 212 The flow is barotropic and the surface moves up and down with gravity as the restoring force. 213 The phase speed of such waves is high (some hundreds of metres per second) so that 214 the time step would have to be very short if they were present in the model. 215 The latter strategy filters out these waves since the rigid lid approximation implies $\eta=0$, 216 $i.e.$ the sea surface is the surface $z=0$. 217 This well known approximation increases the surface wave speed to infinity and 218 modifies certain other longwave dynamics ($e.g.$ barotropic Rossby or planetary waves). 219 The rigid-lid hypothesis is an obsolescent feature in modern OGCMs. 220 It has been available until the release 3.1 of \NEMO, and it has been removed in release 3.2 and followings. 221 Only the free surface formulation is now described in the this document (see the next sub-section). 210 222 211 223 % ------------------------------------------------------------------------------------------------------------- … … 215 227 \label{subsec:PE_free_surface} 216 228 217 In the free surface formulation, a variable $\eta$, the sea-surface height, is introduced218 which describes the shape of the air-sea interface. This variable is solution of a 219 prognostic equation which is established by forming the vertical average of the kinematic 220 surface condition (\autoref{eq:PE_w_bbc}):229 In the free surface formulation, a variable $\eta$, the sea-surface height, 230 is introduced which describes the shape of the air-sea interface. 231 This variable is solution of a prognostic equation which is established by forming the vertical average of 232 the kinematic surface condition (\autoref{eq:PE_w_bbc}): 221 233 \begin{equation} \label{eq:PE_ssh} 222 234 \frac{\partial \eta }{\partial t}=-D+P-E … … 226 238 and using (\autoref{eq:PE_hydrostatic}) the surface pressure is given by: $p_s = \rho \, g \, \eta$. 227 239 228 Allowing the air-sea interface to move introduces the external gravity waves (EGWs) 229 a s a class of solution of the primitive equations. These waves are barotropic because230 of hydrostatic assumption, and their phase speed is quite high. Their time scale is 231 short with respect to the other processes described by the primitive equations.232 233 Two choices can be made regarding the implementation of the free surface in the model, 240 Allowing the air-sea interface to move introduces the external gravity waves (EGWs) as 241 a class of solution of the primitive equations. 242 These waves are barotropic because of hydrostatic assumption, and their phase speed is quite high. 243 Their time scale is short with respect to the other processes described by the primitive equations. 244 245 Two choices can be made regarding the implementation of the free surface in the model, 234 246 depending on the physical processes of interest. 235 247 236 $\bullet$ If one is interested in EGWs, in particular the tides and their interaction 237 with the baroclinic structure of the ocean (internal waves) possibly in shallow seas, 238 then a non linear free surface is the most appropriate. This means that no 239 approximation is made in (\autoref{eq:PE_ssh}) and that the variation of the ocean 240 volume is fully taken into account. Note that in order to study the fast time scales 241 associated with EGWs it is necessary to minimize time filtering effects (use an 242 explicit time scheme with very small time step, or a split-explicit scheme with 243 reasonably small time step, see \autoref{subsec:DYN_spg_exp} or \autoref{subsec:DYN_spg_ts}. 244 245 $\bullet$ If one is not interested in EGW but rather sees them as high frequency 246 noise, it is possible to apply an explicit filter to slow down the fastest waves while 247 not altering the slow barotropic Rossby waves. If further, an approximative conservation 248 of heat and salt contents is sufficient for the problem solved, then it is 249 sufficient to solve a linearized version of (\autoref{eq:PE_ssh}), which still allows 250 to take into account freshwater fluxes applied at the ocean surface \citep{Roullet_Madec_JGR00}. 248 $\bullet$ If one is interested in EGWs, in particular the tides and their interaction with 249 the baroclinic structure of the ocean (internal waves) possibly in shallow seas, 250 then a non linear free surface is the most appropriate. 251 This means that no approximation is made in (\autoref{eq:PE_ssh}) and that 252 the variation of the ocean volume is fully taken into account. 253 Note that in order to study the fast time scales associated with EGWs it is necessary to 254 minimize time filtering effects 255 (use an explicit time scheme with very small time step, or a split-explicit scheme with reasonably small time step, 256 see \autoref{subsec:DYN_spg_exp} or \autoref{subsec:DYN_spg_ts}). 257 258 $\bullet$ If one is not interested in EGW but rather sees them as high frequency noise, 259 it is possible to apply an explicit filter to slow down the fastest waves while 260 not altering the slow barotropic Rossby waves. 261 If further, an approximative conservation of heat and salt contents is sufficient for the problem solved, 262 then it is sufficient to solve a linearized version of (\autoref{eq:PE_ssh}), 263 which still allows to take into account freshwater fluxes applied at the ocean surface \citep{Roullet_Madec_JGR00}. 251 264 Nevertheless, with the linearization, an exact conservation of heat and salt contents is lost. 252 265 253 The filtering of EGWs in models with a free surface is usually a matter of discretisation 254 of the temporal derivatives, using a split-explicit method \citep{Killworth_al_JPO91, Zhang_Endoh_JGR92} 255 or the implicit scheme \citep{Dukowicz1994} or the addition of a filtering force in the momentum equation 256 \citep{Roullet_Madec_JGR00}. With the present release, \NEMO offers the choice between 257 an explicit free surface (see \autoref{subsec:DYN_spg_exp}) or a split-explicit scheme strongly 258 inspired the one proposed by \citet{Shchepetkin_McWilliams_OM05} (see \autoref{subsec:DYN_spg_ts}). 266 The filtering of EGWs in models with a free surface is usually a matter of discretisation of 267 the temporal derivatives, 268 using a split-explicit method \citep{Killworth_al_JPO91, Zhang_Endoh_JGR92} or 269 the implicit scheme \citep{Dukowicz1994} or 270 the addition of a filtering force in the momentum equation \citep{Roullet_Madec_JGR00}. 271 With the present release, \NEMO offers the choice between 272 an explicit free surface (see \autoref{subsec:DYN_spg_exp}) or 273 a split-explicit scheme strongly inspired the one proposed by \citet{Shchepetkin_McWilliams_OM05} 274 (see \autoref{subsec:DYN_spg_ts}). 259 275 260 276 %\newpage … … 274 290 \label{subsec:PE_tensorial} 275 291 276 In many ocean circulation problems, the flow field has regions of enhanced dynamics 277 ($i.e.$ surface layers, western boundary currents, equatorial currents, or ocean fronts). 278 The representation of such dynamical processes can be improved by specifically increasing 279 the model resolution in these regions. As well, it may be convenient to use a lateral 280 boundary-following coordinate system to better represent coastal dynamics. Moreover, 281 the common geographical coordinate system has a singular point at the North Pole that 282 cannot be easily treated in a global model without filtering. A solution consists of introducing 283 an appropriate coordinate transformation that shifts the singular point onto land 284 \citep{Madec_Imbard_CD96, Murray_JCP96}. As a consequence, it is important to solve the primitive 285 equations in various curvilinear coordinate systems. An efficient way of introducing an 286 appropriate coordinate transform can be found when using a tensorial formalism. 287 This formalism is suited to any multidimensional curvilinear coordinate system. Ocean 288 modellers mainly use three-dimensional orthogonal grids on the sphere (spherical earth 289 approximation), with preservation of the local vertical. Here we give the simplified equations 290 for this particular case. The general case is detailed by \citet{Eiseman1980} in their survey 291 of the conservation laws of fluid dynamics. 292 293 Let (\textit{i},\textit{j},\textit{k}) be a set of orthogonal curvilinear coordinates on the sphere 294 associated with the positively oriented orthogonal set of unit vectors (\textbf{i},\textbf{j},\textbf{k}) 295 linked to the earth such that \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are 296 two vectors orthogonal to \textbf{k}, $i.e.$ along geopotential surfaces (\autoref{fig:referential}). 297 Let $(\lambda,\varphi,z)$ be the geographical coordinate system in which a position is defined 298 by the latitude $\varphi(i,j)$, the longitude $\lambda(i,j)$ and the distance from the centre of 299 the earth $a+z(k)$ where $a$ is the earth's radius and $z$ the altitude above a reference sea 300 level (\autoref{fig:referential}). The local deformation of the curvilinear coordinate system is 301 given by $e_1$, $e_2$ and $e_3$, the three scale factors: 292 In many ocean circulation problems, the flow field has regions of enhanced dynamics 293 ($i.e.$ surface layers, western boundary currents, equatorial currents, or ocean fronts). 294 The representation of such dynamical processes can be improved by 295 specifically increasing the model resolution in these regions. 296 As well, it may be convenient to use a lateral boundary-following coordinate system to 297 better represent coastal dynamics. 298 Moreover, the common geographical coordinate system has a singular point at the North Pole that 299 cannot be easily treated in a global model without filtering. 300 A solution consists of introducing an appropriate coordinate transformation that 301 shifts the singular point onto land \citep{Madec_Imbard_CD96, Murray_JCP96}. 302 As a consequence, it is important to solve the primitive equations in various curvilinear coordinate systems. 303 An efficient way of introducing an appropriate coordinate transform can be found when using a tensorial formalism. 304 This formalism is suited to any multidimensional curvilinear coordinate system. 305 Ocean modellers mainly use three-dimensional orthogonal grids on the sphere (spherical earth approximation), 306 with preservation of the local vertical. Here we give the simplified equations for this particular case. 307 The general case is detailed by \citet{Eiseman1980} in their survey of the conservation laws of fluid dynamics. 308 309 Let (\textit{i},\textit{j},\textit{k}) be a set of orthogonal curvilinear coordinates on 310 the sphere associated with the positively oriented orthogonal set of unit vectors 311 (\textbf{i},\textbf{j},\textbf{k}) linked to the earth such that 312 \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are two vectors orthogonal to \textbf{k}, 313 $i.e.$ along geopotential surfaces (\autoref{fig:referential}). 314 Let $(\lambda,\varphi,z)$ be the geographical coordinate system in which a position is defined by 315 the latitude $\varphi(i,j)$, the longitude $\lambda(i,j)$ and 316 the distance from the centre of the earth $a+z(k)$ where $a$ is the earth's radius and 317 $z$ the altitude above a reference sea level (\autoref{fig:referential}). 318 The local deformation of the curvilinear coordinate system is given by $e_1$, $e_2$ and $e_3$, 319 the three scale factors: 302 320 \begin{equation} \label{eq:scale_factors} 303 \begin{aligned} 304 e_1 &=\left( {a+z} \right)\;\left[ {\left( {\frac{\partial \lambda 305 }{\partial i}\cos \varphi } \right)^2+\left( {\frac{\partial \varphi 306 }{\partial i}} \right)^2} \right]^{1/2} \\ 307 e_2 &=\left( {a+z} \right)\;\left[ {\left( {\frac{\partial \lambda 308 }{\partial j}\cos \varphi } \right)^2+\left( {\frac{\partial \varphi 309 }{\partial j}} \right)^2} \right]^{1/2} \\ 310 e_3 &=\left( {\frac{\partial z}{\partial k}} \right) \\ 311 \end{aligned} 312 \end{equation} 313 321 \begin{aligned} 322 e_1 &=\left( {a+z} \right)\;\left[ {\left( {\frac{\partial \lambda}{\partial i}\cos \varphi } \right)^2 323 +\left( {\frac{\partial \varphi }{\partial i}} \right)^2} \right]^{1/2} \\ 324 e_2 &=\left( {a+z} \right)\;\left[ {\left( {\frac{\partial \lambda }{\partial j}\cos \varphi } \right)^2+ 325 \left( {\frac{\partial \varphi }{\partial j}} \right)^2} \right]^{1/2} \\ 326 e_3 &=\left( {\frac{\partial z}{\partial k}} \right) \\ 327 \end{aligned} 328 \end{equation} 329 330 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 331 \begin{figure}[!tb] 332 \begin{center} 333 \includegraphics[width=0.60\textwidth]{Fig_I_earth_referential} 334 \caption{ \protect\label{fig:referential} 335 the geographical coordinate system $(\lambda,\varphi,z)$ and the curvilinear 336 coordinate system (\textbf{i},\textbf{j},\textbf{k}). } 337 \end{center} 338 \end{figure} 314 339 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 315 \begin{figure}[!tb] \begin{center} 316 \includegraphics[width=0.60\textwidth]{Fig_I_earth_referential} 317 \caption{ \protect\label{fig:referential} 318 the geographical coordinate system $(\lambda,\varphi,z)$ and the curvilinear 319 coordinate system (\textbf{i},\textbf{j},\textbf{k}). } 320 \end{center} \end{figure} 321 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 322 323 Since the ocean depth is far smaller than the earth's radius, $a+z$, can be replaced by 324 $a$ in (\autoref{eq:scale_factors}) (thin-shell approximation). The resulting horizontal scale 325 factors $e_1$, $e_2$ are independent of $k$ while the vertical scale factor is a single 326 function of $k$ as \textbf{k} is parallel to \textbf{z}. The scalar and vector operators that 327 appear in the primitive equations (\autoref{eq:PE_dyn} to \autoref{eq:PE_eos}) can 328 be written in the tensorial form, invariant in any orthogonal horizontal curvilinear coordinate 329 system transformation: 340 341 Since the ocean depth is far smaller than the earth's radius, $a+z$, can be replaced by $a$ in 342 (\autoref{eq:scale_factors}) (thin-shell approximation). 343 The resulting horizontal scale factors $e_1$, $e_2$ are independent of $k$ while 344 the vertical scale factor is a single function of $k$ as \textbf{k} is parallel to \textbf{z}. 345 The scalar and vector operators that appear in the primitive equations 346 (\autoref{eq:PE_dyn} to \autoref{eq:PE_eos}) can be written in the tensorial form, 347 invariant in any orthogonal horizontal curvilinear coordinate system transformation: 330 348 \begin{subequations} \label{eq:PE_discrete_operators} 331 349 \begin{equation} \label{eq:PE_grad} … … 369 387 \label{subsec:PE_zco_Eq} 370 388 371 In order to express the Primitive Equations in tensorial formalism, it is necessary to compute 372 the horizontal component of the non-linear and viscous terms of the equation using 373 \autoref{eq:PE_grad}) to \autoref{eq:PE_lap_vector}. 374 Let us set $\vect U=(u,v,w)={\vect{U}}_h +w\;\vect{k}$, the velocity in the $(i,j,k)$ coordinate 375 system and define the relative vorticity $\zeta$ and the divergence of the horizontal velocity 376 field $\chi$, by: 389 In order to express the Primitive Equations in tensorial formalism, 390 it is necessary to compute the horizontal component of the non-linear and viscous terms of the equation using 391 \autoref{eq:PE_grad}) to \autoref{eq:PE_lap_vector}. 392 Let us set $\vect U=(u,v,w)={\vect{U}}_h +w\;\vect{k}$, the velocity in the $(i,j,k)$ coordinate system and 393 define the relative vorticity $\zeta$ and the divergence of the horizontal velocity field $\chi$, by: 377 394 \begin{equation} \label{eq:PE_curl_Uh} 378 395 \zeta =\frac{1}{e_1 e_2 }\left[ {\frac{\partial \left( {e_2 \,v} … … 386 403 \end{equation} 387 404 388 Using the fact that the horizontal scale factors $e_1$ and $e_2$ are independent of $k$ 389 and that $e_3$ is a function of the single variable $k$, the nonlinear term of 390 \autoref{eq:PE_dyn} can be transformed as follows:405 Using the fact that the horizontal scale factors $e_1$ and $e_2$ are independent of $k$ and that 406 $e_3$ is a function of the single variable $k$, 407 the nonlinear term of \autoref{eq:PE_dyn} can be transformed as follows: 391 408 \begin{flalign*} 392 409 &\left[ {\left( { \nabla \times {\rm {\bf U}} } \right) \times {\rm {\bf U}} … … 426 443 \end{flalign*} 427 444 428 The last term of the right hand side is obviously zero, and thus the nonlinear term of 445 The last term of the right hand side is obviously zero, and thus the nonlinear term of 429 446 \autoref{eq:PE_dyn} is written in the $(i,j,k)$ coordinate system: 430 447 \begin{equation} \label{eq:PE_vector_form} … … 437 454 \end{equation} 438 455 439 This is the so-called \textit{vector invariant form} of the momentum advection term. 440 For some purposes, it can be advantageous to write this term in the so-called flux form, 441 $i.e.$ to write it as the divergence of fluxes. For example, the first component of442 \autoref{eq:PE_vector_form} (the $i$-component) is transformed as follows:456 This is the so-called \textit{vector invariant form} of the momentum advection term. 457 For some purposes, it can be advantageous to write this term in the so-called flux form, 458 $i.e.$ to write it as the divergence of fluxes. 459 For example, the first component of \autoref{eq:PE_vector_form} (the $i$-component) is transformed as follows: 443 460 \begin{flalign*} 444 461 &{ \begin{array}{*{20}l} … … 525 542 \end{multline} 526 543 527 The flux form has two terms, the first one is expressed as the divergence of momentum528 fluxes (hence the flux form name given to this formulation) and the second one is due to 529 the curvilinear nature of the coordinate system used. The latter is called the \emph{metric} 530 term and can be viewed as a modification of the Coriolis parameter:544 The flux form has two terms, 545 the first one is expressed as the divergence of momentum fluxes (hence the flux form name given to this formulation) 546 and the second one is due to the curvilinear nature of the coordinate system used. 547 The latter is called the \emph{metric} term and can be viewed as a modification of the Coriolis parameter: 531 548 \begin{equation} \label{eq:PE_cor+metric} 532 549 f \to f + \frac{1}{e_1\;e_2} \left( v \frac{\partial e_2}{\partial i} … … 534 551 \end{equation} 535 552 536 Note that in the case of geographical coordinate, $i.e.$ when $(i,j) \to (\lambda ,\varphi )$537 and $(e_1 ,e_2) \to (a \,\cos \varphi ,a)$, we recover the commonly used modification of 538 the Coriolis parameter $f \to f+(u/a) \tan \varphi$.553 Note that in the case of geographical coordinate, 554 $i.e.$ when $(i,j) \to (\lambda ,\varphi )$ and $(e_1 ,e_2) \to (a \,\cos \varphi ,a)$, 555 we recover the commonly used modification of the Coriolis parameter $f \to f+(u/a) \tan \varphi$. 539 556 540 557 541 558 $\ $\newline % force a new ligne 542 559 543 To sum up, the curvilinear $z$-coordinate equations solved by the ocean model can be 544 written inthe following tensorial formalism:560 To sum up, the curvilinear $z$-coordinate equations solved by the ocean model can be written in 561 the following tensorial formalism: 545 562 546 563 \vspace{+10pt} … … 594 611 \end{multline} 595 612 \end{subequations} 596 where $\zeta$, the relative vorticity, is given by \autoref{eq:PE_curl_Uh} and $p_s $,597 the surface pressure, is given by:613 where $\zeta$, the relative vorticity, is given by \autoref{eq:PE_curl_Uh} and 614 $p_s $, the surface pressure, is given by: 598 615 \begin{equation} \label{eq:PE_spg} 599 616 p_s = \rho \,g \,\eta 600 617 \end{equation} 601 with $\eta$ is solution of \autoref{eq:PE_ssh} 618 with $\eta$ is solution of \autoref{eq:PE_ssh}. 602 619 603 620 The vertical velocity and the hydrostatic pressure are diagnosed from the following equations: … … 628 645 \end{equation} 629 646 630 The expression of \textbf{D}$^{U}$, $D^{S}$ and $D^{T}$ depends on the subgrid scale 631 parameterisation used. It will be defined in \autoref{eq:PE_zdf}. The nature and formulation of 632 ${\rm {\bf F}}^{\rm {\bf U}}$, $F^T$ and $F^S$, the surface forcing terms, are discussed 633 in \autoref{chap:SBC}.647 The expression of \textbf{D}$^{U}$, $D^{S}$ and $D^{T}$ depends on the subgrid scale parameterisation used. 648 It will be defined in \autoref{eq:PE_zdf}. 649 The nature and formulation of ${\rm {\bf F}}^{\rm {\bf U}}$, $F^T$ and $F^S$, the surface forcing terms, 650 are discussed in \autoref{chap:SBC}. 634 651 635 652 … … 642 659 \label{sec:PE_gco} 643 660 644 The ocean domain presents a huge diversity of situation in the vertical. First the ocean surface is a time dependent surface (moving surface). Second the ocean floor depends on the geographical position, varying from more than 6,000 meters in abyssal trenches to zero at the coast. Last but not least, the ocean stratification exerts a strong barrier to vertical motions and mixing. 645 Therefore, in order to represent the ocean with respect to the first point a space and time dependent vertical coordinate that follows the variation of the sea surface height $e.g.$ an $z$*-coordinate; for the second point, a space variation to fit the change of bottom topography $e.g.$ a terrain-following or $\sigma$-coordinate; and for the third point, one will be tempted to use a space and time dependent coordinate that follows the isopycnal surfaces, $e.g.$ an isopycnic coordinate. 646 647 In order to satisfy two or more constrains one can even be tempted to mixed these coordinate systems, as in HYCOM (mixture of $z$-coordinate at the surface, isopycnic coordinate in the ocean interior and $\sigma$ at the ocean bottom) \citep{Chassignet_al_JPO03} or OPA (mixture of $z$-coordinate in vicinity the surface and steep topography areas and $\sigma$-coordinate elsewhere) \citep{Madec_al_JPO96} among others. 648 649 In fact one is totally free to choose any space and time vertical coordinate by introducing an arbitrary vertical coordinate : 661 The ocean domain presents a huge diversity of situation in the vertical. 662 First the ocean surface is a time dependent surface (moving surface). 663 Second the ocean floor depends on the geographical position, 664 varying from more than 6,000 meters in abyssal trenches to zero at the coast. 665 Last but not least, the ocean stratification exerts a strong barrier to vertical motions and mixing. 666 Therefore, in order to represent the ocean with respect to 667 the first point a space and time dependent vertical coordinate that follows the variation of the sea surface height 668 $e.g.$ an $z$*-coordinate; 669 for the second point, a space variation to fit the change of bottom topography 670 $e.g.$ a terrain-following or $\sigma$-coordinate; 671 and for the third point, one will be tempted to use a space and time dependent coordinate that 672 follows the isopycnal surfaces, $e.g.$ an isopycnic coordinate. 673 674 In order to satisfy two or more constrains one can even be tempted to mixed these coordinate systems, as in 675 HYCOM (mixture of $z$-coordinate at the surface, isopycnic coordinate in the ocean interior and $\sigma$ at 676 the ocean bottom) \citep{Chassignet_al_JPO03} or 677 OPA (mixture of $z$-coordinate in vicinity the surface and steep topography areas and $\sigma$-coordinate elsewhere) 678 \citep{Madec_al_JPO96} among others. 679 680 In fact one is totally free to choose any space and time vertical coordinate by 681 introducing an arbitrary vertical coordinate : 650 682 \begin{equation} \label{eq:PE_s} 651 683 s=s(i,j,k,t) 652 684 \end{equation} 653 with the restriction that the above equation gives a single-valued monotonic relationship between $s$ and $k$, when $i$, $j$ and $t$ are held fixed. \autoref{eq:PE_s} is a transformation from the $(i,j,k,t)$ coordinate system with independent variables into the $(i,j,s,t)$ generalised coordinate system with $s$ depending on the other three variables through \autoref{eq:PE_s}. 654 This so-called \textit{generalised vertical coordinate} \citep{Kasahara_MWR74} is in fact an Arbitrary Lagrangian--Eulerian (ALE) coordinate. Indeed, choosing an expression for $s$ is an arbitrary choice that determines which part of the vertical velocity (defined from a fixed referential) will cross the levels (Eulerian part) and which part will be used to move them (Lagrangian part). 655 The coordinate is also sometime referenced as an adaptive coordinate \citep{Hofmeister_al_OM09}, since the coordinate system is adapted in the course of the simulation. Its most often used implementation is via an ALE algorithm, in which a pure lagrangian step is followed by regridding and remapping steps, the later step implicitly embedding the vertical advection \citep{Hirt_al_JCP74, Chassignet_al_JPO03, White_al_JCP09}. Here we follow the \citep{Kasahara_MWR74} strategy : a regridding step (an update of the vertical coordinate) followed by an eulerian step with an explicit computation of vertical advection relative to the moving s-surfaces. 685 with the restriction that the above equation gives a single-valued monotonic relationship between $s$ and $k$, 686 when $i$, $j$ and $t$ are held fixed. 687 \autoref{eq:PE_s} is a transformation from the $(i,j,k,t)$ coordinate system with independent variables into 688 the $(i,j,s,t)$ generalised coordinate system with $s$ depending on the other three variables through 689 \autoref{eq:PE_s}. 690 This so-called \textit{generalised vertical coordinate} \citep{Kasahara_MWR74} is in fact 691 an Arbitrary Lagrangian--Eulerian (ALE) coordinate. 692 Indeed, choosing an expression for $s$ is an arbitrary choice that determines 693 which part of the vertical velocity (defined from a fixed referential) will cross the levels (Eulerian part) and 694 which part will be used to move them (Lagrangian part). 695 The coordinate is also sometime referenced as an adaptive coordinate \citep{Hofmeister_al_OM09}, 696 since the coordinate system is adapted in the course of the simulation. 697 Its most often used implementation is via an ALE algorithm, 698 in which a pure lagrangian step is followed by regridding and remapping steps, 699 the later step implicitly embedding the vertical advection 700 \citep{Hirt_al_JCP74, Chassignet_al_JPO03, White_al_JCP09}. 701 Here we follow the \citep{Kasahara_MWR74} strategy: 702 a regridding step (an update of the vertical coordinate) followed by an eulerian step with 703 an explicit computation of vertical advection relative to the moving s-surfaces. 656 704 657 705 %\gmcomment{ … … 659 707 %A key point here is that the $s$-coordinate depends on $(i,j)$ ==> horizontal pressure gradient... 660 708 661 the generalized vertical coordinates used in ocean modelling are not orthogonal, 662 which contrasts with many other applications in mathematical physics. 663 Hence, it is useful to keep in mind the following properties that may seem 664 odd on initial encounter. 665 666 The horizontal velocity in ocean models measures motions in the horizontal plane, 667 perpendicular to the local gravitational field. That is, horizontal velocity is mathematically 668 the same regardless the vertical coordinate, be it geopotential, isopycnal, pressure, 669 or terrain following. The key motivation for maintaining the same horizontal velocity 670 component is that the hydrostatic and geostrophic balances are dominant in the large-scale ocean. 671 Use of an alternative quasi-horizontal velocity, for example one oriented parallel 672 to the generalized surface, would lead to unacceptable numerical errors. 673 Correspondingly, the vertical direction is anti-parallel to the gravitational force in all 674 of the coordinate systems. We do not choose the alternative of a quasi-vertical 675 direction oriented normal to the surface of a constant generalized vertical coordinate. 676 677 It is the method used to measure transport across the generalized vertical coordinate 678 surfaces which differs between the vertical coordinate choices. That is, computation 679 of the dia-surface velocity component represents the fundamental distinction between 680 the various coordinates. In some models, such as geopotential, pressure, and 681 terrain following, this transport is typically diagnosed from volume or mass conservation. 682 In other models, such as isopycnal layered models, this transport is prescribed based 683 on assumptions about the physical processes producing a flux across the layer interfaces. 684 685 686 In this section we first establish the PE in the generalised vertical $s$-coordinate, 709 the generalized vertical coordinates used in ocean modelling are not orthogonal, 710 which contrasts with many other applications in mathematical physics. 711 Hence, it is useful to keep in mind the following properties that may seem odd on initial encounter. 712 713 The horizontal velocity in ocean models measures motions in the horizontal plane, 714 perpendicular to the local gravitational field. 715 That is, horizontal velocity is mathematically the same regardless the vertical coordinate, be it geopotential, 716 isopycnal, pressure, or terrain following. 717 The key motivation for maintaining the same horizontal velocity component is that 718 the hydrostatic and geostrophic balances are dominant in the large-scale ocean. 719 Use of an alternative quasi-horizontal velocity, for example one oriented parallel to the generalized surface, 720 would lead to unacceptable numerical errors. 721 Correspondingly, the vertical direction is anti-parallel to the gravitational force in 722 all of the coordinate systems. 723 We do not choose the alternative of a quasi-vertical direction oriented normal to 724 the surface of a constant generalized vertical coordinate. 725 726 It is the method used to measure transport across the generalized vertical coordinate surfaces which differs between 727 the vertical coordinate choices. 728 That is, computation of the dia-surface velocity component represents the fundamental distinction between 729 the various coordinates. 730 In some models, such as geopotential, pressure, and terrain following, this transport is typically diagnosed from 731 volume or mass conservation. 732 In other models, such as isopycnal layered models, this transport is prescribed based on assumptions about 733 the physical processes producing a flux across the layer interfaces. 734 735 736 In this section we first establish the PE in the generalised vertical $s$-coordinate, 687 737 then we discuss the particular cases available in \NEMO, namely $z$, $z$*, $s$, and $\tilde z$. 688 738 %} … … 693 743 \subsection{\textit{S-}coordinate formulation} 694 744 695 Starting from the set of equations established in \autoref{sec:PE_zco} for the special case $k=z$ 696 and thus $e_3=1$, we introduce an arbitrary vertical coordinate $s=s(i,j,k,t)$, which includes 697 $z$-, \textit{z*}- and $\sigma-$coordinates as special cases ($s=z$, $s=\textit{z*}$, and 698 $s=\sigma=z/H$ or $=z/\left(H+\eta \right)$, resp.). A formal derivation of the transformed 699 equations is given in \autoref{apdx:A}. Let us define the vertical scale factor by 700 $e_3=\partial_s z$ ($e_3$ is now a function of $(i,j,k,t)$ ), and the slopes in the 701 (\textbf{i},\textbf{j}) directions between $s-$ and $z-$surfaces by:745 Starting from the set of equations established in \autoref{sec:PE_zco} for the special case $k=z$ and thus $e_3=1$, 746 we introduce an arbitrary vertical coordinate $s=s(i,j,k,t)$, 747 which includes $z$-, \textit{z*}- and $\sigma-$coordinates as special cases 748 ($s=z$, $s=\textit{z*}$, and $s=\sigma=z/H$ or $=z/\left(H+\eta \right)$, resp.). 749 A formal derivation of the transformed equations is given in \autoref{apdx:A}. 750 Let us define the vertical scale factor by $e_3=\partial_s z$ ($e_3$ is now a function of $(i,j,k,t)$ ), 751 and the slopes in the (\textbf{i},\textbf{j}) directions between $s-$ and $z-$surfaces by: 702 752 \begin{equation} \label{eq:PE_sco_slope} 703 753 \sigma _1 =\frac{1}{e_1 }\;\left. {\frac{\partial z}{\partial i}} \right|_s … … 711 761 \end{equation} 712 762 713 The equations solved by the ocean model \autoref{eq:PE} in $s-$coordinate can be written as follows (see \autoref{sec:A_momentum}): 763 The equations solved by the ocean model \autoref{eq:PE} in $s-$coordinate can be written as follows 764 (see \autoref{sec:A_momentum}): 714 765 715 766 \vspace{0.5cm} … … 763 814 \end{multline} 764 815 765 where the relative vorticity, \textit{$\zeta $}, the surface pressure gradient, and the hydrostatic766 pressure have the same expressions as in $z$-coordinates although they do not represent 767 exactly the same quantities. $\omega$ is provided by the continuity equation 768 (see \autoref{apdx:A}):816 where the relative vorticity, \textit{$\zeta $}, the surface pressure gradient, 817 and the hydrostatic pressure have the same expressions as in $z$-coordinates although 818 they do not represent exactly the same quantities. 819 $\omega$ is provided by the continuity equation (see \autoref{apdx:A}): 769 820 \begin{equation} \label{eq:PE_sco_continuity} 770 821 \frac{\partial e_3}{\partial t} + e_3 \; \chi + \frac{\partial \omega }{\partial s} = 0 … … 791 842 \end{multline} 792 843 793 The equation of state has the same expression as in $z$-coordinate, and similar expressions794 a re used for mixing and forcing terms.844 The equation of state has the same expression as in $z$-coordinate, 845 and similar expressions are used for mixing and forcing terms. 795 846 796 847 \gmcomment{ … … 808 859 809 860 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 810 \begin{figure}[!b] \begin{center} 811 \includegraphics[width=1.0\textwidth]{Fig_z_zstar} 812 \caption{ \protect\label{fig:z_zstar} 813 (a) $z$-coordinate in linear free-surface case ; 814 (b) $z-$coordinate in non-linear free surface case ; 815 (c) re-scaled height coordinate (become popular as the \textit{z*-}coordinate 816 \citep{Adcroft_Campin_OM04} ).} 817 \end{center} \end{figure} 861 \begin{figure}[!b] 862 \begin{center} 863 \includegraphics[width=1.0\textwidth]{Fig_z_zstar} 864 \caption{ \protect\label{fig:z_zstar} 865 (a) $z$-coordinate in linear free-surface case ; 866 (b) $z-$coordinate in non-linear free surface case ; 867 (c) re-scaled height coordinate (become popular as the \textit{z*-}coordinate 868 \citep{Adcroft_Campin_OM04} ). 869 } 870 \end{center} 871 \end{figure} 818 872 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 819 873 820 874 821 In that case, the free surface equation is nonlinear, and the variations of volume are fully 822 taken into account. These coordinates systems is presented in a report \citep{Levier2007} 823 available on the \NEMO web site. 875 In that case, the free surface equation is nonlinear, and the variations of volume are fully taken into account. 876 These coordinates systems is presented in a report \citep{Levier2007} available on the \NEMO web site. 824 877 825 878 %\gmcomment{ 826 The \textit{z*} coordinate approach is an unapproximated, non-linear free surface implementation 827 which allows one to deal with large amplitude free-surface 828 variations relative to the vertical resolution \citep{Adcroft_Campin_OM04}. In 829 the \textit{z*} formulation, the variation of the column thickness due to sea-surface 830 undulations is not concentrated in the surface level, as in the $z$-coordinate formulation, 831 but is equally distributed over the full water column. Thus vertical 832 levels naturally follow sea-surface variations, with a linear attenuation with 833 depth, as illustrated by figure fig.1c . Note that with a flat bottom, such as in 834 fig.1c, the bottom-following $z$ coordinate and \textit{z*} are equivalent. 835 The definition and modified oceanic equations for the rescaled vertical coordinate 836 \textit{z*}, including the treatment of fresh-water flux at the surface, are 837 detailed in Adcroft and Campin (2004). The major points are summarized 838 here. The position ( \textit{z*}) and vertical discretization (\textit{z*}) are expressed as: 879 The \textit{z*} coordinate approach is an unapproximated, non-linear free surface implementation which allows one to 880 deal with large amplitude free-surface variations relative to the vertical resolution \citep{Adcroft_Campin_OM04}. 881 In the \textit{z*} formulation, 882 the variation of the column thickness due to sea-surface undulations is not concentrated in the surface level, 883 as in the $z$-coordinate formulation, but is equally distributed over the full water column. 884 Thus vertical levels naturally follow sea-surface variations, with a linear attenuation with depth, 885 as illustrated by figure fig.1c. 886 Note that with a flat bottom, such as in fig.1c, the bottom-following $z$ coordinate and \textit{z*} are equivalent. 887 The definition and modified oceanic equations for the rescaled vertical coordinate \textit{z*}, 888 including the treatment of fresh-water flux at the surface, are detailed in Adcroft and Campin (2004). 889 The major points are summarized here. 890 The position ( \textit{z*}) and vertical discretization (\textit{z*}) are expressed as: 839 891 \begin{equation} \label{eq:z-star} 840 892 H + \textit{z*} = (H + z) / r \quad \text{and} \ \delta \textit{z*} = \delta z / r \quad \text{with} \ r = \frac{H+\eta} {H} 841 893 \end{equation} 842 Since the vertical displacement of the free surface is incorporated in the vertical 843 coordinate \textit{z*}, the upper and lower boundaries are at fixed \textit{z*} position, 844 $\textit{z*} = 0$ and $\textit{z*} = -H$ respectively. Also the divergence of the flow field845 is no longer zero as shown by the continuity equation:894 Since the vertical displacement of the free surface is incorporated in the vertical coordinate \textit{z*}, 895 the upper and lower boundaries are at fixed \textit{z*} position, 896 $\textit{z*} = 0$ and $\textit{z*} = -H$ respectively. 897 Also the divergence of the flow field is no longer zero as shown by the continuity equation: 846 898 \begin{equation*} 847 899 \frac{\partial r}{\partial t} = \nabla_{\textit{z*}} \cdot \left( r \; \rm{\bf U}_h \right) … … 853 905 % from MOM4p1 documentation 854 906 855 To overcome problems with vanishing surface and/or bottom cells, we consider the 856 zstar coordinate 907 To overcome problems with vanishing surface and/or bottom cells, we consider the zstar coordinate 857 908 \begin{equation} \label{eq:PE_} 858 909 z^\star = H \left( \frac{z-\eta}{H+\eta} \right) 859 910 \end{equation} 860 911 861 This coordinate is closely related to the "eta" coordinate used in many atmospheric 862 models (see Black (1994) for a review of eta coordinate atmospheric models). It 863 was originally used in ocean models by Stacey et al. (1995) for studies of tides 864 next to shelves, and it has been recently promoted by Adcroft and Campin (2004) 865 for global climate modelling. 866 867 The surfaces of constant $z^\star$ are quasi-horizontal. Indeed, the $z^\star$ coordinate reduces to $z$ when $\eta$ is zero. In general, when noting the large differences between 868 undulations of the bottom topography versus undulations in the surface height, it 869 is clear that surfaces constant $z^\star$ are very similar to the depth surfaces. These properties greatly reduce difficulties of computing the horizontal pressure gradient relative to terrain following sigma models discussed in \autoref{subsec:PE_sco}. 870 Additionally, since $z^\star$ when $\eta = 0$, no flow is spontaneously generated in an 871 unforced ocean starting from rest, regardless the bottom topography. This behaviour is in contrast to the case with "s"-models, where pressure gradient errors in 872 the presence of nontrivial topographic variations can generate nontrivial spontaneous flow from a resting state, depending on the sophistication of the pressure 873 gradient solver. The quasi-horizontal nature of the coordinate surfaces also facilitates the implementation of neutral physics parameterizations in $z^\star$ models using 874 the same techniques as in $z$-models (see Chapters 13-16 of \cite{Griffies_Bk04}) for a 875 discussion of neutral physics in $z$-models, as well as \autoref{sec:LDF_slp} 876 in this document for treatment in \NEMO). 877 878 The range over which $z^\star$ varies is time independent $-H \leq z^\star \leq 0$. Hence, all 879 cells remain nonvanishing, so long as the surface height maintains $\eta > ?H$. This 880 is a minor constraint relative to that encountered on the surface height when using 881 $s = z$ or $s = z - \eta$. 882 883 Because $z^\star$ has a time independent range, all grid cells have static increments 884 ds, and the sum of the ver tical increments yields the time independent ocean 885 depth %·k ds = H. 886 The $z^\star$ coordinate is therefore invisible to undulations of the 887 free surface, since it moves along with the free surface. This proper ty means that 888 no spurious ver tical transpor t is induced across surfaces of constant $z^\star$ by the 889 motion of external gravity waves. Such spurious transpor t can be a problem in 890 z-models, especially those with tidal forcing. Quite generally, the time independent 891 range for the $z^\star$ coordinate is a very convenient proper ty that allows for a nearly 892 arbitrary ver tical resolution even in the presence of large amplitude fluctuations of 893 the surface height, again so long as $\eta > -H$. 912 This coordinate is closely related to the "eta" coordinate used in many atmospheric models 913 (see Black (1994) for a review of eta coordinate atmospheric models). 914 It was originally used in ocean models by Stacey et al. (1995) for studies of tides next to shelves, 915 and it has been recently promoted by Adcroft and Campin (2004) for global climate modelling. 916 917 The surfaces of constant $z^\star$ are quasi-horizontal. 918 Indeed, the $z^\star$ coordinate reduces to $z$ when $\eta$ is zero. 919 In general, when noting the large differences between 920 undulations of the bottom topography versus undulations in the surface height, 921 it is clear that surfaces constant $z^\star$ are very similar to the depth surfaces. 922 These properties greatly reduce difficulties of computing the horizontal pressure gradient relative to 923 terrain following sigma models discussed in \autoref{subsec:PE_sco}. 924 Additionally, since $z^\star$ when $\eta = 0$, 925 no flow is spontaneously generated in an unforced ocean starting from rest, regardless the bottom topography. 926 This behaviour is in contrast to the case with "s"-models, where pressure gradient errors in the presence of 927 nontrivial topographic variations can generate nontrivial spontaneous flow from a resting state, 928 depending on the sophistication of the pressure gradient solver. 929 The quasi-horizontal nature of the coordinate surfaces also facilitates the implementation of 930 neutral physics parameterizations in $z^\star$ models using the same techniques as in $z$-models 931 (see Chapters 13-16 of \cite{Griffies_Bk04}) for a discussion of neutral physics in $z$-models, 932 as well as \autoref{sec:LDF_slp} in this document for treatment in \NEMO). 933 934 The range over which $z^\star$ varies is time independent $-H \leq z^\star \leq 0$. 935 Hence, all cells remain nonvanishing, so long as the surface height maintains $\eta > ?H$. 936 This is a minor constraint relative to that encountered on the surface height when using $s = z$ or $s = z - \eta$. 937 938 Because $z^\star$ has a time independent range, all grid cells have static increments ds, 939 and the sum of the ver tical increments yields the time independent ocean depth. %·k ds = H. 940 The $z^\star$ coordinate is therefore invisible to undulations of the free surface, 941 since it moves along with the free surface. 942 This proper ty means that no spurious vertical transport is induced across surfaces of constant $z^\star$ by 943 the motion of external gravity waves. 944 Such spurious transpor t can be a problem in z-models, especially those with tidal forcing. 945 Quite generally, the time independent range for the $z^\star$ coordinate is a very convenient property that 946 allows for a nearly arbitrary ver tical resolution even in the presence of large amplitude fluctuations of 947 the surface height, again so long as $\eta > -H$. 894 948 895 949 %end MOM doc %%% … … 909 963 \subsubsection{Introduction} 910 964 911 Several important aspects of the ocean circulation are influenced by bottom topography. 912 Of course, the most important is that bottom topography determines deep ocean sub-basins, 913 barriers, sills and channels that strongly constrain the path of water masses, but more subtle 914 effects exist. For example, the topographic $\beta$-effect is usually larger than the planetary 915 one along continental slopes. Topographic Rossby waves can be excited and can interact 916 with the mean current. In the $z-$coordinate system presented in the previous section 917 (\autoref{sec:PE_zco}), $z-$surfaces are geopotential surfaces. The bottom topography is 918 discretised by steps. This often leads to a misrepresentation of a gradually sloping bottom 919 and to large localized depth gradients associated with large localized vertical velocities. 920 The response to such a velocity field often leads to numerical dispersion effects. 921 One solution to strongly reduce this error is to use a partial step representation of bottom 922 topography instead of a full step one \cite{Pacanowski_Gnanadesikan_MWR98}. 923 Another solution is to introduce a terrain-following coordinate system (hereafter $s-$coordinate) 924 925 The $s$-coordinate avoids the discretisation error in the depth field since the layers of 926 computation are gradually adjusted with depth to the ocean bottom. Relatively small 927 topographic features as well as gentle, large-scale slopes of the sea floor in the deep 928 ocean, which would be ignored in typical $z$-model applications with the largest grid 929 spacing at greatest depths, can easily be represented (with relatively low vertical resolution). 930 A terrain-following model (hereafter $s-$model) also facilitates the modelling of the 931 boundary layer flows over a large depth range, which in the framework of the $z$-model 932 would require high vertical resolution over the whole depth range. Moreover, with a 933 $s$-coordinate it is possible, at least in principle, to have the bottom and the sea surface 934 as the only boundaries of the domain (nomore lateral boundary condition to specify). 935 Nevertheless, a $s$-coordinate also has its drawbacks. Perfectly adapted to a 936 homogeneous ocean, it has strong limitations as soon as stratification is introduced. 937 The main two problems come from the truncation error in the horizontal pressure 938 gradient and a possibly increased diapycnal diffusion. The horizontal pressure force 939 in $s$-coordinate consists of two terms (see \autoref{apdx:A}), 965 Several important aspects of the ocean circulation are influenced by bottom topography. 966 Of course, the most important is that bottom topography determines deep ocean sub-basins, barriers, sills and 967 channels that strongly constrain the path of water masses, but more subtle effects exist. 968 For example, the topographic $\beta$-effect is usually larger than the planetary one along continental slopes. 969 Topographic Rossby waves can be excited and can interact with the mean current. 970 In the $z-$coordinate system presented in the previous section (\autoref{sec:PE_zco}), 971 $z-$surfaces are geopotential surfaces. 972 The bottom topography is discretised by steps. 973 This often leads to a misrepresentation of a gradually sloping bottom and to 974 large localized depth gradients associated with large localized vertical velocities. 975 The response to such a velocity field often leads to numerical dispersion effects. 976 One solution to strongly reduce this error is to use a partial step representation of bottom topography instead of 977 a full step one \cite{Pacanowski_Gnanadesikan_MWR98}. 978 Another solution is to introduce a terrain-following coordinate system (hereafter $s-$coordinate). 979 980 The $s$-coordinate avoids the discretisation error in the depth field since the layers of 981 computation are gradually adjusted with depth to the ocean bottom. 982 Relatively small topographic features as well as gentle, large-scale slopes of the sea floor in the deep ocean, 983 which would be ignored in typical $z$-model applications with the largest grid spacing at greatest depths, 984 can easily be represented (with relatively low vertical resolution). 985 A terrain-following model (hereafter $s-$model) also facilitates the modelling of the boundary layer flows over 986 a large depth range, which in the framework of the $z$-model would require high vertical resolution over 987 the whole depth range. 988 Moreover, with a $s$-coordinate it is possible, at least in principle, to have the bottom and the sea surface as 989 the only boundaries of the domain (no more lateral boundary condition to specify). 990 Nevertheless, a $s$-coordinate also has its drawbacks. Perfectly adapted to a homogeneous ocean, 991 it has strong limitations as soon as stratification is introduced. 992 The main two problems come from the truncation error in the horizontal pressure gradient and 993 a possibly increased diapycnal diffusion. 994 The horizontal pressure force in $s$-coordinate consists of two terms (see \autoref{apdx:A}), 940 995 941 996 \begin{equation} \label{eq:PE_p_sco} … … 944 999 \end{equation} 945 1000 946 The second term in \autoref{eq:PE_p_sco} depends on the tilt of the coordinate surface 947 and introduces a truncation error that is not present in a $z$-model. In the special case 948 of a $\sigma-$coordinate (i.e. a depth-normalised coordinate system $\sigma = z/H$), 949 \citet{Haney1991} and \citet{Beckmann1993} have given estimates of the magnitude 950 of this truncation error. It depends on topographic slope, stratification, horizontal and 951 vertical resolution, the equation of state, and the finite difference scheme. This error 952 limits the possible topographic slopes that a model can handle at a given horizontal 953 and vertical resolution. This is a severe restriction for large-scale applications using 954 realistic bottom topography. The large-scale slopes require high horizontal resolution, 955 and the computational cost becomes prohibitive. This problem can be at least partially 956 overcome by mixing $s$-coordinate and step-like representation of bottom topography \citep{Gerdes1993a,Gerdes1993b,Madec_al_JPO96}. However, the definition of the model 957 domain vertical coordinate becomes then a non-trivial thing for a realistic bottom 958 topography: a envelope topography is defined in $s$-coordinate on which a full or 959 partial step bottom topography is then applied in order to adjust the model depth to 960 the observed one (see \autoref{sec:DOM_zgr}. 961 962 For numerical reasons a minimum of diffusion is required along the coordinate surfaces 963 of any finite difference model. It causes spurious diapycnal mixing when coordinate 964 surfaces do not coincide with isoneutral surfaces. This is the case for a $z$-model as 965 well as for a $s$-model. However, density varies more strongly on $s-$surfaces than 966 on horizontal surfaces in regions of large topographic slopes, implying larger diapycnal 967 diffusion in a $s$-model than in a $z$-model. Whereas such a diapycnal diffusion in a 968 $z$-model tends to weaken horizontal density (pressure) gradients and thus the horizontal 969 circulation, it usually reinforces these gradients in a $s$-model, creating spurious circulation. 970 For example, imagine an isolated bump of topography in an ocean at rest with a horizontally 971 uniform stratification. Spurious diffusion along $s$-surfaces will induce a bump of isoneutral 972 surfaces over the topography, and thus will generate there a baroclinic eddy. In contrast, 973 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 974 ($i.e.$ the main thermocline) \citep{Madec_al_JPO96}. An alternate solution consists of rotating 975 the lateral diffusive tensor to geopotential or to isoneutral surfaces (see \autoref{subsec:PE_ldf}. 976 Unfortunately, the slope of isoneutral surfaces relative to the $s$-surfaces can very large, 1001 The second term in \autoref{eq:PE_p_sco} depends on the tilt of the coordinate surface and 1002 introduces a truncation error that is not present in a $z$-model. 1003 In the special case of a $\sigma-$coordinate (i.e. a depth-normalised coordinate system $\sigma = z/H$), 1004 \citet{Haney1991} and \citet{Beckmann1993} have given estimates of the magnitude of this truncation error. 1005 It depends on topographic slope, stratification, horizontal and vertical resolution, the equation of state, 1006 and the finite difference scheme. 1007 This error limits the possible topographic slopes that a model can handle at 1008 a given horizontal and vertical resolution. 1009 This is a severe restriction for large-scale applications using realistic bottom topography. 1010 The large-scale slopes require high horizontal resolution, and the computational cost becomes prohibitive. 1011 This problem can be at least partially overcome by mixing $s$-coordinate and 1012 step-like representation of bottom topography \citep{Gerdes1993a,Gerdes1993b,Madec_al_JPO96}. 1013 However, the definition of the model domain vertical coordinate becomes then a non-trivial thing for 1014 a realistic bottom topography: 1015 a envelope topography is defined in $s$-coordinate on which a full or 1016 partial step bottom topography is then applied in order to adjust the model depth to the observed one 1017 (see \autoref{sec:DOM_zgr}. 1018 1019 For numerical reasons a minimum of diffusion is required along the coordinate surfaces of any finite difference model. 1020 It causes spurious diapycnal mixing when coordinate surfaces do not coincide with isoneutral surfaces. 1021 This is the case for a $z$-model as well as for a $s$-model. 1022 However, density varies more strongly on $s-$surfaces than on horizontal surfaces in regions of 1023 large topographic slopes, implying larger diapycnal diffusion in a $s$-model than in a $z$-model. 1024 Whereas such a diapycnal diffusion in a $z$-model tends to weaken horizontal density (pressure) gradients and thus 1025 the horizontal circulation, it usually reinforces these gradients in a $s$-model, creating spurious circulation. 1026 For example, imagine an isolated bump of topography in an ocean at rest with a horizontally uniform stratification. 1027 Spurious diffusion along $s$-surfaces will induce a bump of isoneutral surfaces over the topography, 1028 and thus will generate there a baroclinic eddy. 1029 In contrast, the ocean will stay at rest in a $z$-model. 1030 As for the truncation error, the problem can be reduced by introducing the terrain-following coordinate below 1031 the strongly stratified portion of the water column ($i.e.$ the main thermocline) \citep{Madec_al_JPO96}. 1032 An alternate solution consists of rotating the lateral diffusive tensor to geopotential or to isoneutral surfaces 1033 (see \autoref{subsec:PE_ldf}). 1034 Unfortunately, the slope of isoneutral surfaces relative to the $s$-surfaces can very large, 977 1035 strongly exceeding the stability limit of such a operator when it is discretized (see \autoref{chap:LDF}). 978 1036 979 The $s-$coordinates introduced here \citep{Lott_al_OM90,Madec_al_JPO96} differ mainly in two 980 aspects from similar models: it allows a representation of bottom topography with mixed 981 full or partial step-like/terrain following topography ; It also offers a completely general 982 transformation, $s=s(i,j,z)$ for the vertical coordinate.1037 The $s-$coordinates introduced here \citep{Lott_al_OM90,Madec_al_JPO96} differ mainly in two aspects from 1038 similar models: 1039 it allows a representation of bottom topography with mixed full or partial step-like/terrain following topography; 1040 It also offers a completely general transformation, $s=s(i,j,z)$ for the vertical coordinate. 983 1041 984 1042 … … 991 1049 992 1050 The $\tilde{z}$-coordinate has been developed by \citet{Leclair_Madec_OM11}. 993 It is available in \NEMO since the version 3.4. Nevertheless, it is currently not robust enough 994 to be used in all possible configurations. Its use is therefore not recommended. 1051 It is available in \NEMO since the version 3.4. 1052 Nevertheless, it is currently not robust enough to be used in all possible configurations. 1053 Its use is therefore not recommended. 995 1054 996 1055 … … 1002 1061 \label{sec:PE_zdf_ldf} 1003 1062 1004 The primitive equations describe the behaviour of a geophysical fluid at 1005 space and time scales larger than a few kilometres in the horizontal, a few 1006 meters in the vertical and a few minutes. They are usually solved at larger 1007 scales: the specified grid spacing and time step of the numerical model. The 1008 effects of smaller scale motions (coming from the advective terms in the 1009 Navier-Stokes equations) must be represented entirely in terms of 1010 large-scale patterns to close the equations. These effects appear in the 1011 equations as the divergence of turbulent fluxes ($i.e.$ fluxes associated with 1012 the mean correlation of small scale perturbations). Assuming a turbulent 1013 closure hypothesis is equivalent to choose a formulation for these fluxes. 1014 It is usually called the subgrid scale physics. It must be emphasized that 1015 this is the weakest part of the primitive equations, but also one of the 1016 most important for long-term simulations as small scale processes \textit{in fine} 1017 balance the surface input of kinetic energy and heat. 1018 1019 The control exerted by gravity on the flow induces a strong anisotropy 1020 between the lateral and vertical motions. Therefore subgrid-scale physics 1021 \textbf{D}$^{\vect{U}}$, $D^{S}$ and $D^{T}$ in \autoref{eq:PE_dyn}, 1022 \autoref{eq:PE_tra_T} and \autoref{eq:PE_tra_S} are divided into a lateral part 1023 \textbf{D}$^{l \vect{U}}$, $D^{lS}$ and $D^{lT}$ and a vertical part 1024 \textbf{D}$^{vU}$, $D^{vS}$ and $D^{vT}$. The formulation of these terms 1025 and their underlying physics are briefly discussed in the next two subsections. 1063 The primitive equations describe the behaviour of a geophysical fluid at space and time scales larger than 1064 a few kilometres in the horizontal, a few meters in the vertical and a few minutes. 1065 They are usually solved at larger scales: the specified grid spacing and time step of the numerical model. 1066 The effects of smaller scale motions (coming from the advective terms in the Navier-Stokes equations) must be represented entirely in terms of large-scale patterns to close the equations. 1067 These effects appear in the equations as the divergence of turbulent fluxes 1068 ($i.e.$ fluxes associated with the mean correlation of small scale perturbations). 1069 Assuming a turbulent closure hypothesis is equivalent to choose a formulation for these fluxes. 1070 It is usually called the subgrid scale physics. 1071 It must be emphasized that this is the weakest part of the primitive equations, 1072 but also one of the most important for long-term simulations as 1073 small scale processes \textit{in fine} balance the surface input of kinetic energy and heat. 1074 1075 The control exerted by gravity on the flow induces a strong anisotropy between the lateral and vertical motions. 1076 Therefore subgrid-scale physics \textbf{D}$^{\vect{U}}$, $D^{S}$ and $D^{T}$ in 1077 \autoref{eq:PE_dyn}, \autoref{eq:PE_tra_T} and \autoref{eq:PE_tra_S} are divided into 1078 a lateral part \textbf{D}$^{l \vect{U}}$, $D^{lS}$ and $D^{lT}$ and 1079 a vertical part \textbf{D}$^{vU}$, $D^{vS}$ and $D^{vT}$. 1080 The formulation of these terms and their underlying physics are briefly discussed in the next two subsections. 1026 1081 1027 1082 % ------------------------------------------------------------------------------------------------------------- … … 1031 1086 \label{subsec:PE_zdf} 1032 1087 1033 The model resolution is always larger than the scale at which the major 1034 sources of vertical turbulence occur (shear instability, internal wave 1035 breaking...). Turbulent motions are thus never explicitly solved, even 1036 partially, but always parameterized. The vertical turbulent fluxes are 1037 assumed to depend linearly on the gradients of large-scale quantities (for 1038 example, the turbulent heat flux is given by $\overline{T'w'}=-A^{vT} \partial_z \overline T$, 1039 where $A^{vT}$ is an eddy coefficient). This formulation is 1040 analogous to that of molecular diffusion and dissipation. This is quite 1041 clearly a necessary compromise: considering only the molecular viscosity 1042 acting on large scale severely underestimates the role of turbulent 1043 diffusion and dissipation, while an accurate consideration of the details of 1044 turbulent motions is simply impractical. The resulting vertical momentum and 1045 tracer diffusive operators are of second order: 1088 The model resolution is always larger than the scale at which the major sources of vertical turbulence occur 1089 (shear instability, internal wave breaking...). 1090 Turbulent motions are thus never explicitly solved, even partially, but always parameterized. 1091 The vertical turbulent fluxes are assumed to depend linearly on the gradients of large-scale quantities 1092 (for example, the turbulent heat flux is given by $\overline{T'w'}=-A^{vT} \partial_z \overline T$, 1093 where $A^{vT}$ is an eddy coefficient). 1094 This formulation is analogous to that of molecular diffusion and dissipation. 1095 This is quite clearly a necessary compromise: considering only the molecular viscosity acting on 1096 large scale severely underestimates the role of turbulent diffusion and dissipation, 1097 while an accurate consideration of the details of turbulent motions is simply impractical. 1098 The resulting vertical momentum and tracer diffusive operators are of second order: 1046 1099 \begin{equation} \label{eq:PE_zdf} 1047 1100 \begin{split} … … 1052 1105 \end{split} 1053 1106 \end{equation} 1054 where $A^{vm}$ and $A^{vT}$ are the vertical eddy viscosity and diffusivity coefficients, 1055 respectively. At the sea surface and at the bottom, turbulent fluxes of momentum, heat 1056 and salt must be specified (see \autoref{chap:SBC} and \autoref{chap:ZDF} and \autoref{sec:TRA_bbl}). 1057 All the vertical physics is embedded in the specification of the eddy coefficients. 1058 They can be assumed to be either constant, or function of the local fluid properties 1059 ($e.g.$ Richardson number, Brunt-Vais\"{a}l\"{a} frequency...), or computed from a 1060 turbulent closure model. The choices available in \NEMO are discussed in \autoref{chap:ZDF}). 1107 where $A^{vm}$ and $A^{vT}$ are the vertical eddy viscosity and diffusivity coefficients, respectively. 1108 At the sea surface and at the bottom, turbulent fluxes of momentum, heat and salt must be specified 1109 (see \autoref{chap:SBC} and \autoref{chap:ZDF} and \autoref{sec:TRA_bbl}). 1110 All the vertical physics is embedded in the specification of the eddy coefficients. 1111 They can be assumed to be either constant, or function of the local fluid properties 1112 ($e.g.$ Richardson number, Brunt-Vais\"{a}l\"{a} frequency...), 1113 or computed from a turbulent closure model. 1114 The choices available in \NEMO are discussed in \autoref{chap:ZDF}). 1061 1115 1062 1116 % ------------------------------------------------------------------------------------------------------------- … … 1066 1120 \label{subsec:PE_ldf} 1067 1121 1068 Lateral turbulence can be roughly divided into a mesoscale turbulence 1069 associated with eddies (which can be solved explicitly if the resolution is 1070 sufficient since their underlying physics are included in the primitive 1071 equations), and a sub mesoscale turbulence which is never explicitly solved 1072 even partially, but always parameterized. The formulation of lateral eddy 1073 fluxes depends on whether the mesoscale is below or above the grid-spacing 1122 Lateral turbulence can be roughly divided into a mesoscale turbulence associated with eddies 1123 (which can be solved explicitly if the resolution is sufficient since 1124 their underlying physics are included in the primitive equations), 1125 and a sub mesoscale turbulence which is never explicitly solved even partially, but always parameterized. 1126 The formulation of lateral eddy fluxes depends on whether the mesoscale is below or above the grid-spacing 1074 1127 ($i.e.$ the model is eddy-resolving or not). 1075 1128 1076 In non-eddy-resolving configurations, the closure is similar to that used 1077 for the vertical physics. The lateral turbulent fluxes are assumed to depend 1078 linearly on the lateral gradients of large-scale quantities. The resulting 1079 lateral diffusive and dissipative operators are of second order. 1080 Observations show that lateral mixing induced by mesoscale turbulence tends 1081 to be along isopycnal surfaces (or more precisely neutral surfaces \cite{McDougall1987}) 1082 rather than across them. 1083 As the slope of neutral surfaces is small in the ocean, a common 1084 approximation is to assume that the `lateral' direction is the horizontal, 1085 $i.e.$ the lateral mixing is performed along geopotential surfaces. This leads 1086 to a geopotential second order operator for lateral subgrid scale physics. 1087 This assumption can be relaxed: the eddy-induced turbulent fluxes can be 1088 better approached by assuming that they depend linearly on the gradients of 1089 large-scale quantities computed along neutral surfaces. In such a case, 1090 the diffusive operator is an isoneutral second order operator and it has 1091 components in the three space directions. However, both horizontal and 1092 isoneutral operators have no effect on mean ($i.e.$ large scale) potential 1093 energy whereas potential energy is a main source of turbulence (through 1094 baroclinic instabilities). \citet{Gent1990} have proposed a 1095 parameterisation of mesoscale eddy-induced turbulence which associates an 1096 eddy-induced velocity to the isoneutral diffusion. Its mean effect is to 1097 reduce the mean potential energy of the ocean. This leads to a formulation 1098 of lateral subgrid-scale physics made up of an isoneutral second order 1099 operator and an eddy induced advective part. In all these lateral diffusive 1100 formulations, the specification of the lateral eddy coefficients remains the 1101 problematic point as there is no really satisfactory formulation of these 1102 coefficients as a function of large-scale features. 1103 1104 In eddy-resolving configurations, a second order operator can be used, but 1105 usually the more scale selective biharmonic operator is preferred as the 1106 grid-spacing is usually not small enough compared to the scale of the 1107 eddies. The role devoted to the subgrid-scale physics is to dissipate the 1108 energy that cascades toward the grid scale and thus to ensure the stability of 1109 the model while not interfering with the resolved mesoscale activity. Another approach 1110 is becoming more and more popular: instead of specifying explicitly a sub-grid scale 1111 term in the momentum and tracer time evolution equations, one uses a advective 1112 scheme which is diffusive enough to maintain the model stability. It must be emphasised 1113 that then, all the sub-grid scale physics is included in the formulation of the 1114 advection scheme. 1115 1116 All these parameterisations of subgrid scale physics have advantages and 1117 drawbacks. There are not all available in \NEMO. For active tracers (temperature and 1118 salinity) the main ones are: Laplacian and bilaplacian operators acting along 1119 geopotential or iso-neutral surfaces, \citet{Gent1990} parameterisation, 1120 and various slightly diffusive advection schemes. 1121 For momentum, the main ones are: Laplacian and bilaplacian operators acting along 1122 geopotential surfaces, and UBS advection schemes when flux form is chosen for the momentum advection. 1129 In non-eddy-resolving configurations, the closure is similar to that used for the vertical physics. 1130 The lateral turbulent fluxes are assumed to depend linearly on the lateral gradients of large-scale quantities. 1131 The resulting lateral diffusive and dissipative operators are of second order. 1132 Observations show that lateral mixing induced by mesoscale turbulence tends to be along isopycnal surfaces 1133 (or more precisely neutral surfaces \cite{McDougall1987}) rather than across them. 1134 As the slope of neutral surfaces is small in the ocean, a common approximation is to assume that 1135 the `lateral' direction is the horizontal, $i.e.$ the lateral mixing is performed along geopotential surfaces. 1136 This leads to a geopotential second order operator for lateral subgrid scale physics. 1137 This assumption can be relaxed: the eddy-induced turbulent fluxes can be better approached by assuming that 1138 they depend linearly on the gradients of large-scale quantities computed along neutral surfaces. 1139 In such a case, the diffusive operator is an isoneutral second order operator and 1140 it has components in the three space directions. 1141 However, 1142 both horizontal and isoneutral operators have no effect on mean ($i.e.$ large scale) potential energy whereas 1143 potential energy is a main source of turbulence (through baroclinic instabilities). 1144 \citet{Gent1990} have proposed a parameterisation of mesoscale eddy-induced turbulence which 1145 associates an eddy-induced velocity to the isoneutral diffusion. 1146 Its mean effect is to reduce the mean potential energy of the ocean. 1147 This leads to a formulation of lateral subgrid-scale physics made up of an isoneutral second order operator and 1148 an eddy induced advective part. 1149 In all these lateral diffusive formulations, 1150 the specification of the lateral eddy coefficients remains the problematic point as 1151 there is no really satisfactory formulation of these coefficients as a function of large-scale features. 1152 1153 In eddy-resolving configurations, a second order operator can be used, 1154 but usually the more scale selective biharmonic operator is preferred as 1155 the grid-spacing is usually not small enough compared to the scale of the eddies. 1156 The role devoted to the subgrid-scale physics is to dissipate the energy that 1157 cascades toward the grid scale and thus to ensure the stability of the model while 1158 not interfering with the resolved mesoscale activity. 1159 Another approach is becoming more and more popular: 1160 instead of specifying explicitly a sub-grid scale term in the momentum and tracer time evolution equations, 1161 one uses a advective scheme which is diffusive enough to maintain the model stability. 1162 It must be emphasised that then, all the sub-grid scale physics is included in the formulation of 1163 the advection scheme. 1164 1165 All these parameterisations of subgrid scale physics have advantages and drawbacks. 1166 There are not all available in \NEMO. For active tracers (temperature and salinity) the main ones are: 1167 Laplacian and bilaplacian operators acting along geopotential or iso-neutral surfaces, 1168 \citet{Gent1990} parameterisation, and various slightly diffusive advection schemes. 1169 For momentum, the main ones are: Laplacian and bilaplacian operators acting along geopotential surfaces, 1170 and UBS advection schemes when flux form is chosen for the momentum advection. 1123 1171 1124 1172 \subsubsection{Lateral laplacian tracer diffusive operator} … … 1133 1181 \end{array} }} \right) 1134 1182 \end{equation} 1135 where $r_1 \;\mbox{and}\;r_2 $ are the slopes between the surface along 1136 which the diffusive operator acts and the model level ($e. g.$ $z$- or 1137 $s$-surfaces). Note that the formulation \autoref{eq:PE_iso_tensor} is exact for the 1138 rotation between geopotential and $s$-surfaces, while it is only an approximation 1139 for the rotation between isoneutral and $z$- or $s$-surfaces. Indeed, in the latter 1140 case, two assumptions are made to simplify \autoref{eq:PE_iso_tensor} \citep{Cox1987}. 1141 First, the horizontal contribution of the dianeutral mixing is neglected since the ratio 1142 between iso and dia-neutral diffusive coefficients is known to be several orders of 1143 magnitude smaller than unity. Second, the two isoneutral directions of diffusion are 1144 assumed to be independent since the slopes are generally less than $10^{-2}$ in the 1145 ocean (see \autoref{apdx:B}). 1146 1147 For \textit{iso-level} diffusion, $r_1$ and $r_2 $ are zero. $\Re $ reduces to the identity 1148 in the horizontal direction, no rotation is applied. 1149 1150 For \textit{geopotential} diffusion, $r_1$ and $r_2 $ are the slopes between the 1151 geopotential and computational surfaces: they are equal to $\sigma _1$ and $\sigma _2$, 1152 respectively (see \autoref{eq:PE_sco_slope}). 1153 1154 For \textit{isoneutral} diffusion $r_1$ and $r_2$ are the slopes between the isoneutral 1155 and computational surfaces. Therefore, they are different quantities, 1156 but have similar expressions in $z$- and $s$-coordinates. In $z$-coordinates: 1183 where $r_1 \;\mbox{and}\;r_2 $ are the slopes between the surface along which the diffusive operator acts and 1184 the model level ($e. g.$ $z$- or $s$-surfaces). 1185 Note that the formulation \autoref{eq:PE_iso_tensor} is exact for 1186 the rotation between geopotential and $s$-surfaces, 1187 while it is only an approximation for the rotation between isoneutral and $z$- or $s$-surfaces. 1188 Indeed, in the latter case, two assumptions are made to simplify \autoref{eq:PE_iso_tensor} \citep{Cox1987}. 1189 First, the horizontal contribution of the dianeutral mixing is neglected since the ratio between iso and 1190 dia-neutral diffusive coefficients is known to be several orders of magnitude smaller than unity. 1191 Second, the two isoneutral directions of diffusion are assumed to be independent since 1192 the slopes are generally less than $10^{-2}$ in the ocean (see \autoref{apdx:B}). 1193 1194 For \textit{iso-level} diffusion, $r_1$ and $r_2 $ are zero. 1195 $\Re $ reduces to the identity in the horizontal direction, no rotation is applied. 1196 1197 For \textit{geopotential} diffusion, 1198 $r_1$ and $r_2 $ are the slopes between the geopotential and computational surfaces: 1199 they are equal to $\sigma _1$ and $\sigma _2$, respectively (see \autoref{eq:PE_sco_slope}). 1200 1201 For \textit{isoneutral} diffusion $r_1$ and $r_2$ are the slopes between the isoneutral and computational surfaces. 1202 Therefore, they are different quantities, but have similar expressions in $z$- and $s$-coordinates. 1203 In $z$-coordinates: 1157 1204 \begin{equation} \label{eq:PE_iso_slopes} 1158 1205 r_1 =\frac{e_3 }{e_1 } \left( \pd[\rho]{i} \right) \left( \pd[\rho]{k} \right)^{-1} \, \quad … … 1162 1209 1163 1210 \subsubsection{Eddy induced velocity} 1164 When the \textit{eddy induced velocity} parametrisation (eiv) \citep{Gent1990} is used, 1211 When the \textit{eddy induced velocity} parametrisation (eiv) \citep{Gent1990} is used, 1165 1212 an additional tracer advection is introduced in combination with the isoneutral diffusion of tracers: 1166 1213 \begin{equation} \label{eq:PE_iso+eiv} … … 1168 1215 +\nabla \cdot \left( {{\vect{U}}^\ast \,T} \right) 1169 1216 \end{equation} 1170 where ${\vect{U}}^\ast =\left( {u^\ast ,v^\ast ,w^\ast } \right)$ is a non-divergent, 1217 where ${\vect{U}}^\ast =\left( {u^\ast ,v^\ast ,w^\ast } \right)$ is a non-divergent, 1171 1218 eddy-induced transport velocity. This velocity field is defined by: 1172 1219 \begin{equation} \label{eq:PE_eiv} … … 1179 1226 \end{split} 1180 1227 \end{equation} 1181 where $A^{eiv}$ is the eddy induced velocity coefficient (or equivalently the isoneutral1182 thickness diffusivity coefficient), and $\tilde{r}_1$ and $\tilde{r}_2$ are the slopes 1183 between isoneutral and \emph{geopotential} surfaces. Their values are 1184 thus independent of the vertical coordinate, but their expression depends on the coordinate:1228 where $A^{eiv}$ is the eddy induced velocity coefficient 1229 (or equivalently the isoneutral thickness diffusivity coefficient), 1230 and $\tilde{r}_1$ and $\tilde{r}_2$ are the slopes between isoneutral and \emph{geopotential} surfaces. 1231 Their values are thus independent of the vertical coordinate, but their expression depends on the coordinate: 1185 1232 \begin{align} \label{eq:PE_slopes_eiv} 1186 1233 \tilde{r}_n = \begin{cases} … … 1191 1238 \end{align} 1192 1239 1193 The normal component of the eddy induced velocity is zero at all the boundaries. 1194 This can be achieved in a model by tapering either the eddy coefficient or the slopes 1195 to zero in the vicinity of the boundaries. The latter strategy is used in \NEMO (cf. \autoref{chap:LDF}). 1240 The normal component of the eddy induced velocity is zero at all the boundaries. 1241 This can be achieved in a model by tapering either the eddy coefficient or 1242 the slopes to zero in the vicinity of the boundaries. 1243 The latter strategy is used in \NEMO (cf. \autoref{chap:LDF}). 1196 1244 1197 1245 \subsubsection{Lateral bilaplacian tracer diffusive operator} … … 1201 1249 D^{lT}= - \Delta \left( \;\Delta T \right) 1202 1250 \qquad \text{where} \;\; \Delta \bullet = \nabla \left( {\sqrt{B^{lT}\,}\;\Re \;\nabla \bullet} \right) 1203 1204 It is the Laplacian operator given by \autoref{eq:PE_iso_tensor} applied twice with 1251 \end{equation} 1252 It is the Laplacian operator given by \autoref{eq:PE_iso_tensor} applied twice with 1205 1253 the harmonic eddy diffusion coefficient set to the square root of the biharmonic one. 1206 1254 … … 1208 1256 \subsubsection{Lateral Laplacian momentum diffusive operator} 1209 1257 1210 The Laplacian momentum diffusive operator along $z$- or $s$-surfaces is found by 1258 The Laplacian momentum diffusive operator along $z$- or $s$-surfaces is found by 1211 1259 applying \autoref{eq:PE_lap_vector} to the horizontal velocity vector (see \autoref{apdx:B}): 1212 1260 \begin{equation} \label{eq:PE_lapU} … … 1224 1272 \end{equation} 1225 1273 1226 Such a formulation ensures a complete separation between the vorticity and 1227 horizontal divergence fields (see \autoref{apdx:C}). 1228 Unfortunately, it is only available in \textit{iso-level} direction. 1229 When a rotation is required ($i.e.$ geopotential diffusion in $s-$coordinates1230 or isoneutral diffusion in both $z$- and $s$-coordinates), the $u$ and $v-$fields 1231 are considered as independent scalar fields, so that the diffusive operator is given by:1274 Such a formulation ensures a complete separation between the vorticity and horizontal divergence fields 1275 (see \autoref{apdx:C}). 1276 Unfortunately, it is only available in \textit{iso-level} direction. 1277 When a rotation is required 1278 ($i.e.$ geopotential diffusion in $s-$coordinates or isoneutral diffusion in both $z$- and $s$-coordinates), 1279 the $u$ and $v-$fields are considered as independent scalar fields, so that the diffusive operator is given by: 1232 1280 \begin{equation} \label{eq:PE_lapU_iso} 1233 1281 \begin{split} 1234 1282 D_u^{l{\rm {\bf U}}} &= \nabla .\left( {A^{lm} \;\Re \;\nabla u} \right) \\ 1235 1283 D_v^{l{\rm {\bf U}}} &= \nabla .\left( {A^{lm} \;\Re \;\nabla v} \right) 1236 \end{split} 1237 \end{equation} 1238 where $\Re$ is given by \autoref{eq:PE_iso_tensor}. It is the same expression as 1239 those used for diffusive operator on tracers. It must be emphasised that such a 1240 formulation is only exact in a Cartesian coordinate system, $i.e.$ on a $f-$ or 1241 $\beta-$plane, not on the sphere. It is also a very good approximation in vicinity 1242 of the Equator in a geographical coordinate system \citep{Lengaigne_al_JGR03}. 1284 \end{split} 1285 \end{equation} 1286 where $\Re$ is given by \autoref{eq:PE_iso_tensor}. 1287 It is the same expression as those used for diffusive operator on tracers. 1288 It must be emphasised that such a formulation is only exact in a Cartesian coordinate system, 1289 $i.e.$ on a $f-$ or $\beta-$plane, not on the sphere. 1290 It is also a very good approximation in vicinity of the Equator in 1291 a geographical coordinate system \citep{Lengaigne_al_JGR03}. 1243 1292 1244 1293 \subsubsection{lateral bilaplacian momentum diffusive operator} 1245 1294 1246 As for tracers, the bilaplacian order momentum diffusive operator is a 1247 re-entering Laplacian operator with the harmonic eddy diffusion coefficient 1248 set to the square root of the biharmonic one. Nevertheless it is currently 1249 not available in the iso-neutral case. 1295 As for tracers, the bilaplacian order momentum diffusive operator is a re-entering Laplacian operator with 1296 the harmonic eddy diffusion coefficient set to the square root of the biharmonic one. 1297 Nevertheless it is currently not available in the iso-neutral case. 1250 1298 1251 1299 \end{document}
Note: See TracChangeset
for help on using the changeset viewer.