/[MITgcm]/manual/s_phys_pkgs/text/gmredi.tex
ViewVC logotype

Contents of /manual/s_phys_pkgs/text/gmredi.tex

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.7 - (show annotations) (download) (as text)
Tue Apr 20 23:32:59 2004 UTC (21 years, 3 months ago) by edhill
Branch: MAIN
Changes since 1.6: +3 -0 lines
File MIME type: application/x-tex
 o add more CMIREDIR tags for updates to the CMI web site

1 \section{Gent/McWiliams/Redi SGS Eddy parameterization}
2 \begin{rawhtml}
3 <!-- CMIREDIR:gmredi: -->
4 \end{rawhtml}
5
6 There are two parts to the Redi/GM parameterization of geostrophic
7 eddies. The first aims to mix tracer properties along isentropes
8 (neutral surfaces) by means of a diffusion operator oriented along the
9 local isentropic surface (Redi). The second part, adiabatically
10 re-arranges tracers through an advective flux where the advecting flow
11 is a function of slope of the isentropic surfaces (GM).
12
13 The first GCM implementation of the Redi scheme was by Cox 1987 in the
14 GFDL ocean circulation model. The original approach failed to
15 distinguish between isopycnals and surfaces of locally referenced
16 potential density (now called neutral surfaces) which are proper
17 isentropes for the ocean. As will be discussed later, it also appears
18 that the Cox implementation is susceptible to a computational mode.
19 Due to this mode, the Cox scheme requires a background lateral
20 diffusion to be present to conserve the integrity of the model fields.
21
22 The GM parameterization was then added to the GFDL code in the form of
23 a non-divergent bolus velocity. The method defines two
24 stream-functions expressed in terms of the isoneutral slopes subject
25 to the boundary condition of zero value on upper and lower
26 boundaries. The horizontal bolus velocities are then the vertical
27 derivative of these functions. Here in lies a problem highlighted by
28 Griffies et al., 1997: the bolus velocities involve multiple
29 derivatives on the potential density field, which can consequently
30 give rise to noise. Griffies et al. point out that the GM bolus fluxes
31 can be identically written as a skew flux which involves fewer
32 differential operators. Further, combining the skew flux formulation
33 and Redi scheme, substantial cancellations take place to the point
34 that the horizontal fluxes are unmodified from the lateral diffusion
35 parameterization.
36
37 \subsection{Redi scheme: Isopycnal diffusion}
38
39 The Redi scheme diffuses tracers along isopycnals and introduces a
40 term in the tendency (rhs) of such a tracer (here $\tau$) of the form:
41 \begin{equation}
42 \bf{\nabla} \cdot \kappa_\rho \bf{K}_{Redi} \bf{\nabla} \tau
43 \end{equation}
44 where $\kappa_\rho$ is the along isopycnal diffusivity and
45 $\bf{K}_{Redi}$ is a rank 2 tensor that projects the gradient of
46 $\tau$ onto the isopycnal surface. The unapproximated projection tensor is:
47 \begin{equation}
48 \bf{K}_{Redi} = \left(
49 \begin{array}{ccc}
50 1 + S_x& S_x S_y & S_x \\
51 S_x S_y & 1 + S_y & S_y \\
52 S_x & S_y & |S|^2 \\
53 \end{array}
54 \right)
55 \end{equation}
56 Here, $S_x = -\partial_x \sigma / \partial_z \sigma$ and $S_y =
57 -\partial_y \sigma / \partial_z \sigma$ are the components of the
58 isoneutral slope.
59
60 The first point to note is that a typical slope in the ocean interior
61 is small, say of the order $10^{-4}$. A maximum slope might be of
62 order $10^{-2}$ and only exceeds such in unstratified regions where
63 the slope is ill defined. It is therefore justifiable, and
64 customary, to make the small slope approximation, $|S| << 1$. The Redi
65 projection tensor then becomes:
66 \begin{equation}
67 \bf{K}_{Redi} = \left(
68 \begin{array}{ccc}
69 1 & 0 & S_x \\
70 0 & 1 & S_y \\
71 S_x & S_y & |S|^2 \\
72 \end{array}
73 \right)
74 \end{equation}
75
76
77 \subsection{GM parameterization}
78
79 The GM parameterization aims to parameterise the ``advective'' or
80 ``transport'' effect of geostrophic eddies by means of a ``bolus''
81 velocity, $\bf{u}^*$. The divergence of this advective flux is added
82 to the tracer tendency equation (on the rhs):
83 \begin{equation}
84 - \bf{\nabla} \cdot \tau \bf{u}^*
85 \end{equation}
86
87 The bolus velocity is defined as:
88 \begin{eqnarray}
89 u^* & = & - \partial_z F_x \\
90 v^* & = & - \partial_z F_y \\
91 w^* & = & \partial_x F_x + \partial_y F_y
92 \end{eqnarray}
93 where $F_x$ and $F_y$ are stream-functions with boundary conditions
94 $F_x=F_y=0$ on upper and lower boundaries. The virtue of casting the
95 bolus velocity in terms of these stream-functions is that they are
96 automatically non-divergent ($\partial_x u^* + \partial_y v^* = -
97 \partial_{xz} F_x - \partial_{yz} F_y = - \partial_z w^*$). $F_x$ and
98 $F_y$ are specified in terms of the isoneutral slopes $S_x$ and $S_y$:
99 \begin{eqnarray}
100 F_x & = & \kappa_{GM} S_x \\
101 F_y & = & \kappa_{GM} S_y
102 \end{eqnarray}
103 This is the form of the GM parameterization as applied by Donabasaglu,
104 1997, in MOM versions 1 and 2.
105
106 \subsection{Griffies Skew Flux}
107
108 Griffies notes that the discretisation of bolus velocities involves
109 multiple layers of differencing and interpolation that potentially
110 lead to noisy fields and computational modes. He pointed out that the
111 bolus flux can be re-written in terms of a non-divergent flux and a
112 skew-flux:
113 \begin{eqnarray*}
114 \bf{u}^* \tau
115 & = &
116 \left( \begin{array}{c}
117 - \partial_z ( \kappa_{GM} S_x ) \tau \\
118 - \partial_z ( \kappa_{GM} S_y ) \tau \\
119 (\partial_x \kappa_{GM} S_x + \partial_y \kappa_{GM} S_y)\tau
120 \end{array} \right)
121 \\
122 & = &
123 \left( \begin{array}{c}
124 - \partial_z ( \kappa_{GM} S_x \tau) \\
125 - \partial_z ( \kappa_{GM} S_y \tau) \\
126 \partial_x ( \kappa_{GM} S_x \tau) + \partial_y ( \kappa_{GM} S_y) \tau)
127 \end{array} \right)
128 + \left( \begin{array}{c}
129 \kappa_{GM} S_x \partial_z \tau \\
130 \kappa_{GM} S_y \partial_z \tau \\
131 - \kappa_{GM} S_x \partial_x \tau - \kappa_{GM} S_y) \partial_y \tau
132 \end{array} \right)
133 \end{eqnarray*}
134 The first vector is non-divergent and thus has no effect on the tracer
135 field and can be dropped. The remaining flux can be written:
136 \begin{equation}
137 \bf{u}^* \tau = - \kappa_{GM} \bf{K}_{GM} \bf{\nabla} \tau
138 \end{equation}
139 where
140 \begin{equation}
141 \bf{K}_{GM} =
142 \left(
143 \begin{array}{ccc}
144 0 & 0 & -S_x \\
145 0 & 0 & -S_y \\
146 S_x & S_y & 0
147 \end{array}
148 \right)
149 \end{equation}
150 is an anti-symmetric tensor.
151
152 This formulation of the GM parameterization involves fewer derivatives
153 than the original and also involves only terms that already appear in
154 the Redi mixing scheme. Indeed, a somewhat fortunate cancellation
155 becomes apparent when we use the GM parameterization in conjunction
156 with the Redi isoneutral mixing scheme:
157 \begin{equation}
158 \kappa_\rho \bf{K}_{Redi} \bf{\nabla} \tau
159 - u^* \tau =
160 ( \kappa_\rho \bf{K}_{Redi} + \kappa_{GM} \bf{K}_{GM} ) \bf{\nabla} \tau
161 \end{equation}
162 In the instance that $\kappa_{GM} = \kappa_{\rho}$ then
163 \begin{equation}
164 \kappa_\rho \bf{K}_{Redi} + \kappa_{GM} \bf{K}_{GM} =
165 \kappa_\rho
166 \left( \begin{array}{ccc}
167 1 & 0 & 0 \\
168 0 & 1 & 0 \\
169 2 S_x & 2 S_y & |S|^2
170 \end{array}
171 \right)
172 \end{equation}
173 which differs from the variable Laplacian diffusion tensor by only
174 two non-zero elements in the $z$-row.
175
176 \fbox{ \begin{minipage}{4.75in}
177 {\em S/R GMREDI\_CALC\_TENSOR} ({\em pkg/gmredi/gmredi\_calc\_tensor.F})
178
179 $\sigma_x$: {\bf SlopeX} (argument on entry)
180
181 $\sigma_y$: {\bf SlopeY} (argument on entry)
182
183 $\sigma_z$: {\bf SlopeY} (argument)
184
185 $S_x$: {\bf SlopeX} (argument on exit)
186
187 $S_y$: {\bf SlopeY} (argument on exit)
188
189 \end{minipage} }
190
191
192
193 \subsection{Variable $\kappa_{GM}$}
194
195 Visbeck et al., 1996, suggest making the eddy coefficient,
196 $\kappa_{GM}$, a function of the Eady growth rate,
197 $|f|/\sqrt{Ri}$. The formula involves a non-dimensional constant,
198 $\alpha$, and a length-scale $L$:
199 \begin{displaymath}
200 \kappa_{GM} = \alpha L^2 \overline{ \frac{|f|}{\sqrt{Ri}} }^z
201 \end{displaymath}
202 where the Eady growth rate has been depth averaged (indicated by the
203 over-line). A local Richardson number is defined $Ri = N^2 / (\partial
204 u/\partial z)^2$ which, when combined with thermal wind gives:
205 \begin{displaymath}
206 \frac{1}{Ri} = \frac{(\frac{\partial u}{\partial z})^2}{N^2} =
207 \frac{ ( \frac{g}{f \rho_o} | {\bf \nabla} \sigma | )^2 }{N^2} =
208 \frac{ M^4 }{ |f|^2 N^2 }
209 \end{displaymath}
210 where $M^2$ is defined $M^2 = \frac{g}{\rho_o} |{\bf \nabla} \sigma|$.
211 Substituting into the formula for $\kappa_{GM}$ gives:
212 \begin{displaymath}
213 \kappa_{GM} = \alpha L^2 \overline{ \frac{M^2}{N} }^z =
214 \alpha L^2 \overline{ \frac{M^2}{N^2} N }^z =
215 \alpha L^2 \overline{ |S| N }^z
216 \end{displaymath}
217
218
219 \subsection{Tapering and stability}
220
221 Experience with the GFDL model showed that the GM scheme has to be
222 matched to the convective parameterization. This was originally
223 expressed in connection with the introduction of the KPP boundary
224 layer scheme (Large et al., 97) but in fact, as subsequent experience
225 with the MIT model has found, is necessary for any convective
226 parameterization.
227
228 \fbox{ \begin{minipage}{4.75in}
229 {\em S/R GMREDI\_SLOPE\_LIMIT} ({\em
230 pkg/gmredi/gmredi\_slope\_limit.F})
231
232 $\sigma_x, s_x$: {\bf SlopeX} (argument)
233
234 $\sigma_y, s_y$: {\bf SlopeY} (argument)
235
236 $\sigma_z$: {\bf dSigmadRReal} (argument)
237
238 $z_\sigma^{*}$: {\bf dRdSigmaLtd} (argument)
239
240 \end{minipage} }
241
242 \begin{figure}
243 \begin{center}
244 \resizebox{5.0in}{3.0in}{\includegraphics{part6/tapers.eps}}
245 \end{center}
246 \caption{Taper functions used in GKW99 and DM95.}
247 \label{fig:tapers}
248 \end{figure}
249
250 \begin{figure}
251 \begin{center}
252 \resizebox{5.0in}{3.0in}{\includegraphics{part6/effective_slopes.eps}}
253 \end{center}
254 \caption{Effective slope as a function of ``true'' slope using Cox
255 slope clipping, GKW91 limiting and DM95 limiting.}
256 \label{fig:effective_slopes}
257 \end{figure}
258
259
260 \subsubsection{Slope clipping}
261
262 Deep convection sites and the mixed layer are indicated by
263 homogenized, unstable or nearly unstable stratification. The slopes in
264 such regions can be either infinite, very large with a sign reversal
265 or simply very large. From a numerical point of view, large slopes
266 lead to large variations in the tensor elements (implying large bolus
267 flow) and can be numerically unstable. This was first recognized by
268 Cox, 1987, who implemented ``slope clipping'' in the isopycnal mixing
269 tensor. Here, the slope magnitude is simply restricted by an upper
270 limit:
271 \begin{eqnarray}
272 |\nabla \sigma| & = & \sqrt{ \sigma_x^2 + \sigma_y^2 } \\
273 S_{lim} & = & - \frac{|\nabla \sigma|}{ S_{max} }
274 \;\;\;\;\;\;\;\; \mbox{where $S_{max}$ is a parameter} \\
275 \sigma_z^\star & = & \min( \sigma_z , S_{lim} ) \\
276 {[s_x,s_y]} & = & - \frac{ [\sigma_x,\sigma_y] }{\sigma_z^\star}
277 \end{eqnarray}
278 Notice that this algorithm assumes stable stratification through the
279 ``min'' function. In the case where the fluid is well stratified ($\sigma_z < S_{lim}$) then the slopes evaluate to:
280 \begin{equation}
281 {[s_x,s_y]} = - \frac{ [\sigma_x,\sigma_y] }{\sigma_z}
282 \end{equation}
283 while in the limited regions ($\sigma_z > S_{lim}$) the slopes become:
284 \begin{equation}
285 {[s_x,s_y]} = \frac{ [\sigma_x,\sigma_y] }{|\nabla \sigma|/S_{max}}
286 \end{equation}
287 so that the slope magnitude is limited $\sqrt{s_x^2 + s_y^2} =
288 S_{max}$.
289
290 The slope clipping scheme is activated in the model by setting {\bf
291 GM\_tap\-er\_scheme = 'clipping'} in {\em data.gmredi}.
292
293 Even using slope clipping, it is normally the case that the vertical
294 diffusion term (with coefficient $\kappa_\rho{\bf K}_{33} =
295 \kappa_\rho S_{max}^2$) is large and must be time-stepped using an
296 implicit procedure (see section on discretisation and code later).
297 Fig. \ref{fig-mixedlayer} shows the mixed layer depth resulting from
298 a) using the GM scheme with clipping and b) no GM scheme (horizontal
299 diffusion). The classic result of dramatically reduced mixed layers is
300 evident. Indeed, the deep convection sites to just one or two points
301 each and are much shallower than we might prefer. This, it turns out,
302 is due to the over zealous re-stratification due to the bolus transport
303 parameterization. Limiting the slopes also breaks the adiabatic nature
304 of the GM/Redi parameterization, re-introducing diabatic fluxes in
305 regions where the limiting is in effect.
306
307 \subsubsection{Tapering: Gerdes, Koberle and Willebrand, Clim. Dyn. 1991}
308
309 The tapering scheme used in Gerdes et al., 1999, (\cite{gkw:99})
310 addressed two issues with the clipping method: the introduction of
311 large vertical fluxes in addition to convective adjustment fluxes is
312 avoided by tapering the GM/Redi slopes back to zero in
313 low-stratification regions; the adjustment of slopes is replaced by a
314 tapering of the entire GM/Redi tensor. This means the direction of
315 fluxes is unaffected as the amplitude is scaled.
316
317 The scheme inserts a tapering function, $f_1(S)$, in front of the
318 GM/Redi tensor:
319 \begin{equation}
320 f_1(S) = \min \left[ 1, \left( \frac{S_{max}}{|S|}\right)^2 \right]
321 \end{equation}
322 where $S_{max}$ is the maximum slope you want allowed. Where the
323 slopes, $|S|<S_{max}$ then $f_1(S) = 1$ and the tensor is un-tapered
324 but where $|S| \ge S_{max}$ then $f_1(S)$ scales down the tensor so
325 that the effective vertical diffusivity term $\kappa f_1(S) |S|^2 =
326 \kappa S_{max}^2$.
327
328 The GKW tapering scheme is activated in the model by setting {\bf
329 GM\_tap\-er\_scheme = 'gkw91'} in {\em data.gmredi}.
330
331 \subsection{Tapering: Danabasoglu and McWilliams, J. Clim. 1995}
332
333 The tapering scheme used by Danabasoglu and McWilliams, 1995,
334 \cite{dm:95}, followed a similar procedure but used a different
335 tapering function, $f_1(S)$:
336 \begin{equation}
337 f_1(S) = \frac{1}{2} \left( 1+\tanh \left[ \frac{S_c - |S|}{S_d} \right] \right)
338 \end{equation}
339 where $S_c = 0.004$ is a cut-off slope and $S_d=0.001$ is a scale over
340 which the slopes are smoothly tapered. Functionally, the operates in
341 the same way as the GKW91 scheme but has a substantially lower
342 cut-off, turning off the GM/Redi SGS parameterization for weaker
343 slopes.
344
345 The DM tapering scheme is activated in the model by setting {\bf
346 GM\_tap\-er\_scheme = 'dm95'} in {\em data.gmredi}.
347
348 \subsection{Tapering: Large, Danabasoglu and Doney, JPO 1997}
349
350 The tapering used in Large et al., 1997, \cite{ldd:97}, is based on the
351 DM95 tapering scheme, but also tapers the scheme with an additional
352 function of height, $f_2(z)$, so that the GM/Redi SGS fluxes are
353 reduced near the surface:
354 \begin{equation}
355 f_2(S) = \frac{1}{2} \left( 1 + \sin(\pi \frac{z}{D} - \pi/2)\right)
356 \end{equation}
357 where $D = L_\rho |S|$ is a depth-scale and $L_\rho=c/f$ with
358 $c=2$~m~s$^{-1}$. This tapering with height was introduced to fix
359 some spurious interaction with the mixed-layer KPP parameterization.
360
361 The LDD tapering scheme is activated in the model by setting {\bf
362 GM\_tap\-er\_scheme = 'ldd97'} in {\em data.gmredi}.
363
364
365
366
367 \begin{figure}
368 \begin{center}
369 %\includegraphics{mixedlayer-cox.eps}
370 %\includegraphics{mixedlayer-diff.eps}
371 Figure missing.
372 \end{center}
373 \caption{Mixed layer depth using GM parameterization with a) Cox slope
374 clipping and for comparison b) using horizontal constant diffusion.}
375 \label{fig-mixedlayer}
376 \end{figure}
377
378 \subsection{Package Reference}
379 % DO NOT EDIT HERE
380
381
382

  ViewVC Help
Powered by ViewVC 1.1.22