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