meanflowIntro.tex
10.2 KB

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231

```
%
%$Id: meanflowIntro.tex,v 1.4 2005-08-16 14:42:17 hb Exp $
%
\section{The mean flow model \label{sec:meanflowIntro}}
\subsection{Introduction}
This module contains the definitions of the most important
mean flow variables used in geophysical models. In GOTM, these
are
\begin{itemize}
\item the mean horizontal velocity components, $U$ and $V$
\item the mean potential temperature, $\Theta$, (or the mean buoyancy, $B$)
\item the mean salinity, $S$
\end{itemize}
Note that in general a variable $\phi$ describing a turbulent
field can be decomposed into a mean and a fluctuating part. In GOTM,
we use the notation
\begin{equation}
\label{decomposition}
\phi = \mean{\phi} + \phi'
\comma
\end{equation}
where $\mean{\rule{3mm}{0mm}}$ denotes the ensemble mean and the prime
the fluctuating part. In addition, for brevity, we use the following conventions:
\begin{equation}
\label{decompositionConventions}
\begin{array}{rcl}
U &=& \mean{u} \quad \text{for the x-velocity} \\
V &=& \mean{v} \quad \text{for the y-velocity} \\
P &=& \mean{p} \quad \text{for the pressure} \\
\Theta &=& \mean{\theta} \quad \text{for the potential temperature} \\
B &=& \mean{b} \quad \text{for the buoyancy} \\
S &=& \mean{s} \quad \text{for the salinity}
\end{array}
\end{equation}
Note that, if not explicitly mentioned, GOTM uses the units kg, m, s,
K. Further conventions are introduced in the turbulence chapter
\sect{sec:turbulenceIntro}. All operations on these meanflow variables
are executed and coordinated in the {\tt meanflow} module.
\subsubsection{Physics}\label{sec:meanflowIntroPhysics}
Due to the one-dimensional character of GOTM, the state-variables
listed above are assumed to be horizontally homogeneous, depending
only on the vertical $z$-coordinate. As a consequence, all
horizontal gradients have to be taken from observations, or they have
to be estimated, parameterised or neglected.
For example, the surface slopes $\partial_x\zeta$ and
$\partial_y\zeta$ representing the barotropic pressure-gradients may
be determined by means of local observations or results from
three-dimensional numerical models. It is also possible to prescribe a
time series of the near-bed velocity components for reconstructing the
barotropic pressure gradient, see \cite{Burchard99}. The
implementation of these options for the external pressure gradient is
carried out in {\tt extpressure.F90}, described in
\sect{sec:extpressure}. The internal pressure-gradient, which results
from horizontal density gradients, can be prescribed from observations
of horizontal gradients of $\Theta$ and $S$ or from three-dimensional
model results (see {\tt intpressure.F90} in \sect{sec:intpressure}).
These gradients may also be used for horizontally advecting $\Theta$
and $S$ (see \sect{sec:temperature} and \sect{sec:salinity}).
Another option in GOTM for parameterising the advection of $\Theta$
and $S$ is to relax the model results to observations. Evidently, this
raises questions about the physical consistency of the model, but it
might help to provide a more realistic density field for studies of
turbulence dynamics. Nudging is also possible for the horizontal
velocity components. This makes sense in order to initialise inertial
oscillations from observed velocity profiles, see \sect{sec:uequation}
and \sect{sec:vequation}. In the momentum equations, advection and
horizontal diffusion terms are neglected.
In hydrostatic ocean models, the vertical velocity is calculated by
means of the continuity equation, where the horizontal gradients of
$U$ and $V$ are needed. Since these are not available or set to zero,
the assumption of zero vertical velocity would be consistent. In many
applications however, a non-zero vertical velocity is needed in order
to reflect the vertical adiabatic motion of e.g.\ a thermocline. In
GOTM, we have thus included the option of prescribing a vertical
velocity time series at one height level which might be vertically
moving. Vertical velocities at the surface and at the bottom are
prescribed according to the kinematic boundary conditions ($w=0$ at
the bottom and $w=\partial_t\zeta$ at the surface), and between these
locations and the prescribed vertical velocity at a certain height,
linear interpolation is applied, see {\tt updategrid.F90} in
\sect{sec:updategrid}. This vertical velocity is then used for the
vertical advection of all prognostic quantities.
Standard relations according to the law of the wall are used for
deriving bottom boundary conditions for the momentum equations (see
{\tt friction.F90} in \sect{sec:friction}). At the sea surface, they
have to be prescribed or calculated from meteorological observations
with the aid of bulk formulae using the simulated or observed sea
surface temperature (see \sect{sec:airsea}). In {\tt
stratification.F90} described in \sect{sec:stratification}, the
buoyancy $b$ as defined in equation \eq{DefBuoyancy} is calculated by
means of the UNESCO equation of state (\cite{FofonoffMillard83})
or its linearised version. In
special cases, the buoyancy may also be calculated from a simple
transport equation. {\tt stratification.F90} is also used for
calculating the Brunt-V\"ais\"al\"a frequency, $N$.
The turbulent fluxes are calculated by means of various different
turbulence closure models described in great detail in the {\tt
turbulence} module, see \sect{sec:turbulence}. As a simplifying
alternative, mixing can be computed according to the so-called
`convective adjustment' algorithm, see \sect{sec:convective}.
Furthermore, the vertical grid is also defined in the meanflow module
(see {\tt updategrid.F90} in \sect{sec:updategrid}). Choices for the
numerical grid are so-called $\sigma$-coordinates with layers heights
having a fixed portion of the water depth throughout the
simulation. Equidistant and non-equidistant grids are possible.
\subsubsection{Numerics}\label{SectionNumericsMean}
For the spatial discretisation, the water column is divided into $N_i$
layers of not necessarily equal thickness $h_i$,
\begin{equation}
\label{grid}
h_i=(\gamma_i-\gamma_{i-1})D, \qquad i=1,\dots,N_i
\comma
\end{equation}
with nondimensional interfaces $\gamma_i$ with $\gamma_0=-1$,
$\gamma_{i-1}< \gamma_i$ and $\gamma_{N_i}=0$,
see \cite{BurchardPetersen97}.
The discrete values for the mean flow quantities $U$, $V$, $\Theta$,
and $S$ represent interval means and are therefore located at the
centres of the intervals, and the turbulent quantities like $k$, $L$,
$\epsilon$, $\nu_t$, $\nu'_t$, $N$, $P$, $G$, $c_{\mu}$, and
$c_{\mu}'$ are positioned at the interfaces of the intervals (see
\sect{sec:turbulence}). The indexing is such, that the interface
above an interval has the same index as the interval itself. This
means that mean flow quantities range from $i=1,..,N_i$ while
turbulent quantities range from $i=0,..,N_i$ (see \fig{FigGrid}).
\begin{figure}[!h]
\begin{center}
\scalebox{0.5}{\includegraphics{figures/gridvert.eps}}
\caption{Spatial organisation and indexing of the numerical grid.\label{FigGrid}}
\end{center}
\end{figure}
The staggering of the grid allows for a
straight-forward discretisation of the vertical fluxes of momentum
and tracers without averaging. However, for the vertical fluxes of
e.g.\ $k$ and $\epsilon$, averaging of the eddy diffusivities is
necessary. This is only problematic for the fluxes near the surface
and the bottom, where viscosities at the boundaries have to be
considered for the averaging. These can however be derived from the
law of the wall.
\begin{figure}
\begin{center}
\scalebox{0.5}{\includegraphics{figures/gridtime.eps}}
\caption{Temporal organisation and indexing of the numerical grid.
Here, a time stepping slightly more implicit than the
\index{Crank-Nicolson scheme} \cite{CrankNicolson47}
scheme with $\sigma=0.6$ is shown.\label{FigGridTime}}
\end{center}
\end{figure}
The time stepping is equidistant, based on two time levels and not
limited by Courant numbers, because of the absence of advection and an
implicit treatment of vertical diffusion, see \fig{FigGridTime}. In
the following, the discretisation of a simple diffusion equation,
\begin{equation}
\label{simpleDiffusion}
\partder{X}{t} - \partder{}{z} \left( \nu \partder{X}{z} \right) = 0
\comma
\end{equation}
will be illustrated for Neumann-type
boundary conditions
\begin{equation}
\nu \partder{X}{z} = F_s
\qquad \mbox{for } z=\zeta,\qquad
\end{equation}
and
\begin{equation}
\nu \partder{X}{z} = F_b
\qquad \mbox{for } z=-H.\qquad
\end{equation}
The semi-implicit discretisation of \eq{simpleDiffusion}
can then be written as
\begin{equation}\label{sigmafirst}
\displaystyle
\frac{X^{n+1}_{N_i}-X^n_{N_i}}{\Delta t}
-\frac{F_s
-\nu^n_{N_i-1}\frac{X^{n+\sigma}_{N_i}-X^{n+\sigma}_{N_i-1}}{0.5(h^{n+1}_{N_i}+h^{n+1}_{N_i-1})}}{h^{n
+1}_{N_i}}
=
\comma
\end{equation}
\begin{equation}\label{Xdiscrete}
\displaystyle
\frac{X^{n+1}_i-X^n_i}{\Delta t}
-\frac{\nu^n_i\frac{X^{n+\sigma}_{i+1}-X^{n+\sigma}_{i}}{0.5(h^{n+1}_{i+1}+h^{n+1}_i)}
-\nu^n_{i-1}\frac{X^{n+\sigma}_{i}-X^{n+\sigma}_{i-1}}{0.5(h^{n+1}_i+h^{n+1}_{i-1})}}{h^{n
+1}_i}
=0
\comma
\end{equation}
\begin{equation}\label{sigmalast}
\displaystyle
\frac{X^{n+1}_1-X^n_1}{\Delta t}
-\frac{\nu^n_1\frac{X^{n+\sigma}_{2}-X^{n+\sigma}_{1}}{0.5(h^{n+1}_{2}+h^{n+1}_1)}
-F_b}{h^{n+1}_1}
=0
\comma
\end{equation}
for $1<i<N_i$. Here, the semi-implicit time level is defined by
\begin{equation}
X^{n+\sigma}=\sigma X^{n+1}+(1-\sigma)X^n.
\end{equation}
Thus, for $\sigma=0$, a fully explicit, for $\sigma=1$ a fully
implicit, and for $\sigma=0.5$ the \cite{CrankNicolson47}
second-order scheme are obtained. \Fig{FigGridTime} shows an
example for $\sigma=0.6$. It should be noted that often a time
stepping is preferable which is slightly more implicit than the
\cite{CrankNicolson47} scheme in order to obtain
asymptotic stability. The resulting linear system of equations
(\ref{sigmafirst}) -- (\ref{sigmalast}) with tri-diagonal matrix
structure is solved by means of the simplified Gaussian elimination.
With the same strategy, a very similar system of equations can be
derived for variables located at the interfaces of the grid cells,
i.e. variables describing turbulence.
```