- 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_conservation.tex
r9407 r10354 9 9 \minitoc 10 10 11 The continuous equations of motion have many analytic properties. Many 12 quantities (total mass, energy, enstrophy, etc.) are strictly conserved in 13 the inviscid and unforced limit, while ocean physics conserve the total 14 quantities on which they act (momentum, temperature, salinity) but dissipate 15 their total variance (energy, enstrophy, etc.). Unfortunately, the finite 16 difference form of these equations is not guaranteed to retain all these 17 important properties. In constructing the finite differencing schemes, we 18 wish to ensure that certain integral constraints will be maintained. In 19 particular, it is desirable to construct the finite difference equations so 20 that horizontal kinetic energy and/or potential enstrophy of horizontally 21 non-divergent flow, and variance of temperature and salinity will be 22 conserved in the absence of dissipative effects and forcing. \citet{Arakawa1966} 23 has first pointed out the advantage of this approach. He showed that if 24 integral constraints on energy are maintained, the computation will be free 25 of the troublesome "non linear" instability originally pointed out by 26 \citet{Phillips1959}. A consistent formulation of the energetic properties is 27 also extremely important in carrying out long-term numerical simulations for 28 an oceanographic model. Such a formulation avoids systematic errors that 29 accumulate with time \citep{Bryan1997}. 30 31 The general philosophy of OPA which has led to the discrete formulation 32 presented in {\S}II.2 and II.3 is to choose second order non-diffusive 33 scheme for advective terms for both dynamical and tracer equations. At this 34 level of complexity, the resulting schemes are dispersive schemes. 35 Therefore, they require the addition of a diffusive operator to be stable. 36 The alternative is to use diffusive schemes such as upstream or flux 37 corrected schemes. This last option was rejected because we prefer a 38 complete handling of the model diffusion, i.e. of the model physics rather 39 than letting the advective scheme produces its own implicit diffusion 40 without controlling the space and time structure of this implicit diffusion. 41 Note that in some very specific cases as passive tracer studies, the 42 positivity of the advective scheme is required. In that case, and in that 43 case only, the advective scheme used for passive tracer is a flux correction 44 scheme \citep{Marti1992, Levy1996, Levy1998}. 11 The continuous equations of motion have many analytic properties. 12 Many quantities (total mass, energy, enstrophy, etc.) are strictly conserved in the inviscid and unforced limit, 13 while ocean physics conserve the total quantities on which they act (momentum, temperature, salinity) but 14 dissipate their total variance (energy, enstrophy, etc.). 15 Unfortunately, the finite difference form of these equations is not guaranteed to 16 retain all these important properties. 17 In constructing the finite differencing schemes, we wish to ensure that 18 certain integral constraints will be maintained. 19 In particular, it is desirable to construct the finite difference equations so that 20 horizontal kinetic energy and/or potential enstrophy of horizontally non-divergent flow, 21 and variance of temperature and salinity will be conserved in the absence of dissipative effects and forcing. 22 \citet{Arakawa1966} has first pointed out the advantage of this approach. 23 He showed that if integral constraints on energy are maintained, 24 the computation will be free of the troublesome "non linear" instability originally pointed out by 25 \citet{Phillips1959}. 26 A consistent formulation of the energetic properties is also extremely important in carrying out 27 long-term numerical simulations for an oceanographic model. 28 Such a formulation avoids systematic errors that accumulate with time \citep{Bryan1997}. 29 30 The general philosophy of OPA which has led to the discrete formulation presented in {\S}II.2 and II.3 is to 31 choose second order non-diffusive scheme for advective terms for both dynamical and tracer equations. 32 At this level of complexity, the resulting schemes are dispersive schemes. 33 Therefore, they require the addition of a diffusive operator to be stable. 34 The alternative is to use diffusive schemes such as upstream or flux corrected schemes. 35 This last option was rejected because we prefer a complete handling of the model diffusion, 36 i.e. of the model physics rather than letting the advective scheme produces its own implicit diffusion without 37 controlling the space and time structure of this implicit diffusion. 38 Note that in some very specific cases as passive tracer studies, the positivity of the advective scheme is required. 39 In that case, and in that case only, the advective scheme used for passive tracer is a flux correction scheme 40 \citep{Marti1992, Levy1996, Levy1998}. 45 41 46 42 % ------------------------------------------------------------------------------------------------------------- … … 50 46 \label{sec:Invariant_dyn} 51 47 52 The non linear term of the momentum equations has been split into a 53 vorticity term, a gradient of horizontal kinetic energy and a vertical 54 advection term. Three schemes are available for the former (see {\S}~II.2) 55 according to the CPP variable defined (default option\textbf{ 56 }or \textbf{key{\_}vorenergy } or \textbf{key{\_}vorcombined 57 } defined). They differ in their conservative 58 properties (energy or enstrophy conserving scheme). The two latter terms 59 preserve the total kinetic energy: the large scale kinetic energy is also 60 preserved in practice. The remaining non-diffusive terms of the momentum 61 equation (namely the hydrostatic and surface pressure gradient terms) also 62 preserve the total kinetic energy and have no effect on the vorticity of the 63 flow. 48 The non linear term of the momentum equations has been split into a vorticity term, 49 a gradient of horizontal kinetic energy and a vertical advection term. 50 Three schemes are available for the former (see {\S}~II.2) according to the CPP variable defined 51 (default option\textbf{?}or \textbf{key{\_}vorenergy} or \textbf{key{\_}vorcombined} defined). 52 They differ in their conservative properties (energy or enstrophy conserving scheme). 53 The two latter terms preserve the total kinetic energy: 54 the large scale kinetic energy is also preserved in practice. 55 The remaining non-diffusive terms of the momentum equation 56 (namely the hydrostatic and surface pressure gradient terms) also preserve the total kinetic energy and 57 have no effect on the vorticity of the flow. 64 58 65 59 \textbf{* relative, planetary and total vorticity term:} 66 60 67 Let us define as either the relative, planetary and total potential 68 vorticity, i.e. , , and , respectively. The continuous formulation of the 69 vorticity term satisfies following integral constraints: 61 Let us define as either the relative, planetary and total potential vorticity, i.e. ?, ?, and ?, respectively. 62 The continuous formulation of the vorticity term satisfies following integral constraints: 70 63 \begin{equation} \label{eq:vor_vorticity} 71 64 \int_D {{\textbf {k}}\cdot \frac{1}{e_3 }\nabla \times \left( {\varsigma … … 82 75 \int_D {{\textbf{U}}_h \times \left( {\varsigma \;{\textbf{k}}\times {\textbf{U}}_h } \right)\;dv} =0 83 76 \end{equation} 84 where $dv = e_1\, e_2\, e_3\, di\, dj\, dk$ is the volume element. 85 (II.4.1a) means that $\varsigma $ is conserved. (II.4.1b) is obtained by an 86 integration by part. It means that $\varsigma^2$ is conserved for a horizontally 87 non-divergent flow. 88 (II.4.1c) is even satisfied locally since the vorticity term is orthogonal 89 to the horizontal velocity. It means that the vorticity term has no 90 contribution to the evolution of the total kinetic energy. (II.4.1a) is 91 obviously always satisfied, but (II.4.1b) and (II.4.1c) cannot be satisfied 92 simultaneously with a second order scheme. Using the symmetry or 93 anti-symmetry properties of the operators (Eqs II.1.10 and 11), it can be 94 shown that the scheme (II.2.11) satisfies (II.4.1b) but not (II.4.1c), while 95 scheme (II.2.12) satisfies (II.4.1c) but not (II.4.1b) (see appendix C). 96 Note that the enstrophy conserving scheme on total vorticity has been chosen 97 as the standard discrete form of the vorticity term. 77 where $dv = e_1\, e_2\, e_3\, di\, dj\, dk$ is the volume element. 78 (II.4.1a) means that $\varsigma $ is conserved. (II.4.1b) is obtained by an integration by part. 79 It means that $\varsigma^2$ is conserved for a horizontally non-divergent flow. 80 (II.4.1c) is even satisfied locally since the vorticity term is orthogonal to the horizontal velocity. 81 It means that the vorticity term has no contribution to the evolution of the total kinetic energy. 82 (II.4.1a) is obviously always satisfied, but (II.4.1b) and (II.4.1c) cannot be satisfied simultaneously with 83 a second order scheme. 84 Using the symmetry or anti-symmetry properties of the operators (Eqs II.1.10 and 11), 85 it can be shown that the scheme (II.2.11) satisfies (II.4.1b) but not (II.4.1c), 86 while scheme (II.2.12) satisfies (II.4.1c) but not (II.4.1b) (see appendix C). 87 Note that the enstrophy conserving scheme on total vorticity has been chosen as the standard discrete form of 88 the vorticity term. 98 89 99 90 \textbf{* Gradient of kinetic energy / vertical advection} 100 91 101 In continuous formulation, the gradient of horizontal kinetic energy has no 102 contribution to the evolution of the vorticity as the curl of a gradient is 103 zero. This property is satisfied locally with the discrete form of both the 104 gradient and the curl operator we have made (property (II.1.9)~). Another 105 continuous property is that the change of horizontal kinetic energy due to 106 vertical advection is exactly balanced by the change of horizontal kinetic 107 energy due tothe horizontal gradient of horizontal kinetic energy:92 In continuous formulation, the gradient of horizontal kinetic energy has no contribution to the evolution of 93 the vorticity as the curl of a gradient is zero. 94 This property is satisfied locally with the discrete form of both the gradient and the curl operator we have made 95 (property (II.1.9)~). 96 Another continuous property is that the change of horizontal kinetic energy due to 97 vertical advection is exactly balanced by the change of horizontal kinetic energy due to 98 the horizontal gradient of horizontal kinetic energy: 108 99 109 100 \begin{equation} \label{eq:keg_zad} … … 112 103 \end{equation} 113 104 114 Using the discrete form given in {\S}II.2-a and the symmetry or 115 anti-symmetry properties of the mean and difference operators, \autoref{eq:keg_zad} is 116 demonstrated in the Appendix C. The main point here is that satisfying 117 \autoref{eq:keg_zad} links the choice of the discrete forms of the vertical advection 118 and of the horizontal gradient of horizontal kinetic energy. Choosing one 119 imposes the other. The discrete form of the vertical advection given in 120 {\S}II.2-a is a direct consequence of formulating the horizontal kinetic 121 energy as $1/2 \left( \overline{u^2}^i + \overline{v^2}^j \right) $ inthe gradient term.105 Using the discrete form given in {\S}II.2-a and the symmetry or anti-symmetry properties of 106 the mean and difference operators, \autoref{eq:keg_zad} is demonstrated in the Appendix C. 107 The main point here is that satisfying \autoref{eq:keg_zad} links the choice of the discrete forms of 108 the vertical advection and of the horizontal gradient of horizontal kinetic energy. 109 Choosing one imposes the other. 110 The discrete form of the vertical advection given in {\S}II.2-a is a direct consequence of 111 formulating the horizontal kinetic energy as $1/2 \left( \overline{u^2}^i + \overline{v^2}^j \right) $ in 112 the gradient term. 122 113 123 114 \textbf{* hydrostatic pressure gradient term} 124 115 125 In continuous formulation, a pressure gradient has no contribution to the 126 evolution of the vorticity as the curl of a gradient is zero. This 127 properties is satisfied locally with the choice of discretization we have 128 made (property (II.1.9)~). In addition, when the equation of state is linear 129 (i.e. when an advective-diffusive equation for density can be derived from 130 those of temperature and salinity) the change of horizontal kinetic energy 131 due to the work of pressure forces is balanced by the change of potential 132 energy due to buoyancy forces: 116 In continuous formulation, a pressure gradient has no contribution to the evolution of the vorticity as 117 the curl of a gradient is zero. 118 This properties is satisfied locally with the choice of discretization we have made (property (II.1.9)~). 119 In addition, when the equation of state is linear 120 (i.e. when an advective-diffusive equation for density can be derived from those of temperature and salinity) 121 the change of horizontal kinetic energy due to the work of pressure forces is balanced by the change of 122 potential energy due to buoyancy forces: 133 123 134 124 \begin{equation} \label{eq:hpg_pe} … … 136 126 \end{equation} 137 127 138 Using the discrete form given in {\S}~II.2-a and the symmetry or 139 anti-symmetry properties of the mean and difference operators, (II.4.3) is 140 demonstrated in the Appendix C. The main point here is that satisfying 141 (II.4.3) strongly constraints the discrete expression of the depth of 142 $T$-points and of the term added to the pressure gradient in $s-$coordinates: the 143 depth of a $T$-point, $z_T$, is defined as the sum the vertical scale 144 factors at $w$-points starting from the surface. 128 Using the discrete form given in {\S}~II.2-a and the symmetry or anti-symmetry properties of 129 the mean and difference operators, (II.4.3) is demonstrated in the Appendix C. 130 The main point here is that satisfying (II.4.3) strongly constraints the discrete expression of the depth of 131 $T$-points and of the term added to the pressure gradient in $s-$coordinates: the depth of a $T$-point, $z_T$, 132 is defined as the sum the vertical scale factors at $w$-points starting from the surface. 145 133 146 134 \textbf{* surface pressure gradient term} 147 135 148 In continuous formulation, the surface pressure gradient has no contribution 149 to the evolution of vorticity. This properties is trivially satisfied 150 locally as (II.2.3) (the equation verified by $\psi$ has been 151 derived from the discrete formulation of the momentum equations, vertical 152 sum and curl. Nevertheless, the $\psi$-equation is solved numerically by an 153 iterative solver (see {\S}~III.5), thus the property is only satisfied with 154 the accuracy required on the solver. In addition, with the rigid-lid 155 approximation, the change of horizontal kinetic energy due to the work of 136 In continuous formulation, the surface pressure gradient has no contribution to the evolution of vorticity. 137 This properties is trivially satisfied locally as (II.2.3) 138 (the equation verified by $\psi$ has been derived from the discrete formulation of the momentum equations, 139 vertical sum and curl). 140 Nevertheless, the $\psi$-equation is solved numerically by an iterative solver (see {\S}~III.5), 141 thus the property is only satisfied with the accuracy required on the solver. 142 In addition, with the rigid-lid approximation, the change of horizontal kinetic energy due to the work of 156 143 surface pressure forces is exactly zero: 157 144 \begin{equation} \label{eq:spg} … … 159 146 \end{equation} 160 147 161 (II.4.4) is satisfied in discrete form only if the discrete barotropic 162 streamfunction time evolution equation is given by (II.2.3) (see appendix 163 C). This shows that (II.2.3) is the only way to compute the streamfunction, 164 otherwise there is no guarantee that the surface pressure force work 165 vanishes. 148 (II.4.4) is satisfied in discrete form only if 149 the discrete barotropic streamfunction time evolution equation is given by (II.2.3) (see appendix C). 150 This shows that (II.2.3) is the only way to compute the streamfunction, 151 otherwise there is no guarantee that the surface pressure force work vanishes. 166 152 167 153 % ------------------------------------------------------------------------------------------------------------- … … 171 157 \label{sec:Invariant_tra} 172 158 173 In continuous formulation, the advective terms of the tracer equations 174 conserve the tracer content and the quadratic form of the tracer, i.e. 159 In continuous formulation, the advective terms of the tracer equations conserve the tracer content and 160 the quadratic form of the tracer, $i.e.$ 175 161 \begin{equation} \label{eq:tra_tra2} 176 162 \int_D {\nabla .\left( {T\;{\textbf{U}}} \right)\;dv} =0 … … 179 165 \end{equation} 180 166 181 The numerical scheme used ({\S}II.2-b) (equations in flux form, second order 182 centred finite differences) satisfies (II.4.5) (see appendix C). Note that 183 in both continuous and discrete formulations, there is generally no strict 184 conservation of mass, since the equation of state is non linear with respect 185 to $T$ and $S$.In practice, the mass is conserved with a very good accuracy.167 The numerical scheme used ({\S}II.2-b) (equations in flux form, second order centred finite differences) satisfies 168 (II.4.5) (see appendix C). 169 Note that in both continuous and discrete formulations, there is generally no strict conservation of mass, 170 since the equation of state is non linear with respect to $T$ and $S$. 171 In practice, the mass is conserved with a very good accuracy. 186 172 187 173 % ------------------------------------------------------------------------------------------------------------- … … 193 179 \textbf{* lateral momentum diffusion term} 194 180 195 The continuous formulation of the horizontal diffusion of momentum satisfies 196 the following integral constraints~: 181 The continuous formulation of the horizontal diffusion of momentum satisfies the following integral constraints~: 197 182 \begin{equation} \label{eq:dynldf_dyn} 198 183 \int\limits_D {\frac{1}{e_3 }{\rm {\bf k}}\cdot \nabla \times \left[ {\nabla … … 227 212 228 213 229 (II.4.6a) and (II.4.6b) means that the horizontal diffusion of momentum 230 conserve both the potential vorticity and the divergence of the flow, while 231 Eqs (II.4.6c) to (II.4.6e) mean that it dissipates the energy, the enstrophy 232 and the square of the divergence. The two latter properties are only 233 satisfied when the eddy coefficients are horizontally uniform. 234 235 Using (II.1.8) and (II.1.9), and the symmetry or anti-symmetry properties of 236 the mean and difference operators, it is shown that the discrete form of the 237 lateral momentum diffusion given in {\S}II.2-c satisfies all the integral 238 constraints (II.4.6) (see appendix C). In particular, when the eddy 239 coefficients are horizontally uniform, a complete separation of vorticity 240 and horizontal divergence fields is ensured, so that diffusion (dissipation) 241 of vorticity (enstrophy) does not generate horizontal divergence (variance 242 of the horizontal divergence) and \textit{vice versa}. When the vertical curl of the horizontal 243 diffusion of momentum (discrete sense) is taken, the term associated to the 244 horizontal gradient of the divergence is zero locally. When the horizontal 245 divergence of the horizontal diffusion of momentum (discrete sense) is 246 taken, the term associated to the vertical curl of the vorticity is zero 247 locally. The resulting term conserves $\chi$ and dissipates 248 $\chi^2$ when the 249 eddy coefficient is horizontally uniform. 214 (II.4.6a) and (II.4.6b) means that the horizontal diffusion of momentum conserve both the potential vorticity and 215 the divergence of the flow, while Eqs (II.4.6c) to (II.4.6e) mean that it dissipates the energy, the enstrophy and 216 the square of the divergence. 217 The two latter properties are only satisfied when the eddy coefficients are horizontally uniform. 218 219 Using (II.1.8) and (II.1.9), and the symmetry or anti-symmetry properties of the mean and difference operators, 220 it is shown that the discrete form of the lateral momentum diffusion given in 221 {\S}II.2-c satisfies all the integral constraints (II.4.6) (see appendix C). 222 In particular, when the eddy coefficients are horizontally uniform, 223 a complete separation of vorticity and horizontal divergence fields is ensured, 224 so that diffusion (dissipation) of vorticity (enstrophy) does not generate horizontal divergence 225 (variance of the horizontal divergence) and \textit{vice versa}. 226 When the vertical curl of the horizontal diffusion of momentum (discrete sense) is taken, 227 the term associated to the horizontal gradient of the divergence is zero locally. 228 When the horizontal divergence of the horizontal diffusion of momentum (discrete sense) is taken, 229 the term associated to the vertical curl of the vorticity is zero locally. 230 The resulting term conserves $\chi$ and dissipates $\chi^2$ when the eddy coefficient is horizontally uniform. 250 231 251 232 \textbf{* vertical momentum diffusion term} 252 233 253 As for the lateral momentum physics, the continuous form of the vertical 254 diffusion ofmomentum satisfies following integral constraints~:234 As for the lateral momentum physics, the continuous form of the vertical diffusion of 235 momentum satisfies following integral constraints~: 255 236 256 237 conservation of momentum, dissipation of horizontal kinetic energy … … 273 254 \end{aligned} 274 255 \end{equation} 275 conservation of horizontal divergence, dissipation of square of the 276 horizontal divergence 256 conservation of horizontal divergence, dissipation of square of the horizontal divergence 277 257 \begin{equation} \label{eq:dynzdf_div} 278 258 \begin{aligned} … … 286 266 \end{equation} 287 267 288 In discrete form, all these properties are satisfied in $z$-coordinate (see 289 Appendix C). In $s$-coordinates, only first order properties can be 290 demonstrated, i.e. the vertical momentum physics conserve momentum, 291 potential vorticity, and horizontal divergence. 268 In discrete form, all these properties are satisfied in $z$-coordinate (see Appendix C). 269 In $s$-coordinates, only first order properties can be demonstrated, 270 $i.e.$ the vertical momentum physics conserve momentum, potential vorticity, and horizontal divergence. 292 271 293 272 % ------------------------------------------------------------------------------------------------------------- … … 297 276 \label{subsec:Invariant_tra_physics} 298 277 299 The numerical schemes used for tracer subgridscale physics are written in 300 such a way that the heat and salt contents are conserved (equations in flux 301 form, second order centred finite differences). As a form flux is used to 302 compute the temperature and salinity, the quadratic form of these quantities 303 (i.e. their variance) globally tends to diminish. As for the advective term, 304 there is generally no strict conservation of mass even if, in practice, the 305 mass is conserved with a very good accuracy. 278 The numerical schemes used for tracer subgridscale physics are written in such a way that 279 the heat and salt contents are conserved (equations in flux form, second order centred finite differences). 280 As a form flux is used to compute the temperature and salinity, 281 the quadratic form of these quantities (i.e. their variance) globally tends to diminish. 282 As for the advective term, there is generally no strict conservation of mass even if, 283 in practice, the mass is conserved with a very good accuracy. 306 284 307 285 \textbf{* lateral physics: }conservation of tracer, dissipation of tracer … … 315 293 \end{equation} 316 294 317 \textbf{* vertical physics: }conservation of tracer, dissipation of tracer 318 variance, i.e. 295 \textbf{* vertical physics: }conservation of tracer, dissipation of tracer variance, $i.e.$ 319 296 320 297 \begin{equation} \label{eq:trazdf_t_t2} … … 325 302 \end{equation} 326 303 327 Using the symmetry or anti-symmetry properties of the mean and difference 328 operators, it is shown that the discrete form of tracer physics given in 329 {\S}~II.2-c satisfies all the integral constraints (II.4.8) and (II.4.9) 330 except the dissipation of the square of the tracer when non-geopotential 331 diffusion is used (see appendix C). A discrete form of the lateral tracer 332 physics can be derived which satisfies these last properties. Nevertheless, 333 it requires a horizontal averaging of the vertical component of the lateral 334 physics that prevents the use of implicit resolution in the vertical. It has 335 not been implemented. 304 Using the symmetry or anti-symmetry properties of the mean and difference operators, 305 it is shown that the discrete form of tracer physics given in {\S}~II.2-c satisfies all the integral constraints 306 (II.4.8) and (II.4.9) except the dissipation of the square of the tracer when non-geopotential diffusion is used 307 (see appendix C). 308 A discrete form of the lateral tracer physics can be derived which satisfies these last properties. 309 Nevertheless, it requires a horizontal averaging of the vertical component of the lateral physics that 310 prevents the use of implicit resolution in the vertical. 311 It has not been implemented. 336 312 337 313 \end{document}
Note: See TracChangeset
for help on using the changeset viewer.
