Horizontal diffusion

Horizontal diffusion can be very efficient at suppressing small scale noise ($2d$ wavelength) which can quickly contaminate a numerical simulation. In addition, physically-based approaches attempt to represent the cascade of energy owing to turbulent mixing. However, a balance between a minimum suppression of noise and prevention of over-smoothing (or even in the limit, numerical instabilities) must be achieved. There are several options for horizontal diffusion. One can use 2nd, 4th or 6th order diffusion: the second order ($o=2$) diffusion uses a variable deformation-dependent diffusivity which is based on the classic so-called Smagorinsky (1963) model. Hyper (or high order) diffusion ($o$= 4 or 6) is available with either a constant diffusivity (available for 4th and 6th order diffusion) or using the appropriately scaled deformation-dependent diffusivity (available for 4th order diffusion). High-order diffusion is generally regarded as a filter to remove small scale noise which might develop (while leaving horizontal gradients relatively in tact).

The horizontal diffusion of order $o$ for a quantity $\varphi$ can be expressed as

$$D_{h\,\varphi}^o = {\frac{\partial F^{o-1}_{\varphi\,x}}{\partial x}} \,+\, {\frac{\partial F^{o-1}_{\varphi\,y}}{\partial y}}$$ (90)

where the diffusive flux is expressed as

$$F^{o-1}_{\varphi\,x} = \wp_o\,\kappa_\varphi {\frac{\partial^{o-1} \varphi}{\partial x^{o-1}}} \hskip1.375in (\varphi=u,\, v,\, q_{cl},\,q_{ci},\,e )$$ (91)

$$F^{o-1}_{\varphi\,x} = \wp_o\,\kappa_\varphi \left( {\frac{\partial^{o-1} \varphi}{\pa... ...\partial^{o-1} \varphi}{\partial x^{o-1}}}\right) \hskip.2in (\varphi=T,\, q_v)$$ (92)

Note that the tendencies of $T$ and $q_j$ are converted to $\theta_\rho$ and $r_j$ tendencies using the conversion factors detailed in previous sections. The leading term on the RHS of Eq. 91 is defined in Eq. 86 (it is a toggle to determine the sign). The horizontal geopotential gradients in Eq. 93 take into account the height difference (departure) from the reference ($\eta$) surface (i.e. they are used to interpolate the variables to a quasi-horizontal surface). This is done because interpolation of these variables without adjustment (i.e. along coordinate surfaces) can result in errors where the surfaces are steeply sloped. $\Gamma_\varphi$ represents the vertical lapse rate of the variable of interest ($\varphi$) which is defined as

$$\Gamma_T = \omega_\Gamma \, {\frac{\partial T}{\partial z}} - \left(1-\omega_\Gamma\right)\Gamma_0 \hskip.4in \left(\Gamma_d \leq \Gamma_T \leq 0\right)$$ (93)

$$\Gamma_q = \omega_\Gamma \, {\frac{\partial q}{\partial z}} - \left(1-\ome... ..._0{\frac{\partial q_{sat}}{\partial T}} \hskip.5in \left(\Gamma_q \leq 0\right)$$ (94)

where $\Gamma_d$ and $\Gamma_0$ represent the dry adiabatic (-0.010 K m$^{-1}$) and standard (-0.0065 K m$^{-1}$) atmospheric temperature lapse rates, respectively, and $q_{sat}$ represents the saturation specific humidity (kg kg$^{-1}$). The weight factor is

$$\omega_\Gamma = {\frac {1}{1+s_\eta \left(\nabla\varphi_s/\varphi_c \right)}} \hskip.5in \left( 1 \geq \omega_\Gamma \geq 0 \right)$$ (95)

and the gradient of the surface geopotential is defined as

$$\nabla\varphi_s = {\rm max} \left({ \Big\vert {\frac{\partial\var... ...ig\vert ,\, \Big\vert {\frac{\partial\varphi_s}{\partial y}} \Big\vert }\right)$$ (96)

The $\varphi_c$ term in Eq. 96 is a scaling factor defined as 0.001$g$ (where $g$ is the gravitational constant). This formulation is based on the work of F. Giorgi and Nieman. (1993) who used this coefficient to reduce horizontal diffusion in the presence of strong topographical gradients. In ASP, there is also a vertical dependence included which is not present in F. Giorgi and Nieman. (1993), where

$$s_\eta = {\rm min}\left(\eta, \, \eta_b\right) \hskip.5in \left( 0 \leq s_\eta \leq 1 \right)$$ (97)

The coefficient $\eta_b$ is defined as

$$\eta_b={\frac{\pi_m-\pi_t}{p_0-\pi_t}} \hskip.5in \left( 0 \leq \eta_b \leq 1 \right)$$ (98)

Note that $\eta_b$ is computed in the pre-time integration part of the model using a standard pressure profile. The vertical dependence of $s_\eta$ (and therefore $\omega_\Gamma$) causes the weight to increase with increasing altitude (as the coordinate surface begin to flatten out relative to surfaces near the surface), since this reduction factor is mainly needed in the lower atmosphere where the terrain-following coordinate surfaces are steeply sloped in regions of sharp topographic gradients. The result is that the weight $\omega_\Gamma$ is generally unity in upper levels of the atmosphere and in regions with low topographic variability, so that the lapse rates in Eq. 93 are determined by the actual local lapse rates (the general case). But note that in regions of high topographic variability, the lapse rates can vary significantly horizontally, thus introducing anomalous cooling/heating moistening/drying. Some GCM models combat this by using globally averaged lapse rates. Here, we take a somewhat simpler but similar approach by using a standard lapse rate, which all but eliminates the errors discussed above. Some models zero the horizontal diffusion in regions of steep topography (e.g.s F. Giorgi and Nieman. (1993); Zangl (2002)) but one can argue that the high order diffusion schemes act as a high order noise filter, it is not desirable for the diffusion to vanish in regions of large topographic variations. In ASP, both options exist. More details of this parameterization will be given in the section on the vertical discretization. This is a computationally efficient and fairly robust alternative to actually computing true horizontal diffusion in terrain following coordinates.

The diffusivity in Eq. 91 for second order diffusion ($o=2$) of some quantity $\varphi$ is expressed as

$$\kappa_\varphi = \kappa_h \hskip1.00in (\varphi \,=\, u,\,v)$$ (99)

$$= \kappa_h/P_r \hskip0.80in (\varphi \,=\, T,\,q_v,\,q_{cl},\,q_{ci},\,e)$$ (100)

where the turbulent Prandtl number (here defined as 1/3) is used for mass variables. $\kappa_h$ (m$^2$ s$^{-1}$) represents the spatially and temporally varying deformation-dependent (second order) Smagorinsky-type diffusivity (the parameterization is described in the Physics section).

An example of second order diffusion of $\varphi$ from Eq. 91 is

$$D_{h\,\varphi}^2 = {\frac{1}{\mu_d}}\left[ {\frac{\partial}{\part... ...left({ {\kappa_{\varphi}} {\frac{\partial\varphi}{\partial y}} }\right) \right]$$ (101)

For higher order (or hyper) diffusion ($o=4$ or $o=6$), the diffusion coefficient is constant and is defined

$$\kappa_\varphi^o = {\frac{k_h \, d^{o} \, \alpha_\varphi }{2^{o} \, \Delta t}}$$ (102)

The $\kappa_\varphi^o$ in Eq. 103 utilizes a scaled form of the second-order physically-based Smagorinsky coefficient. Hyper diffusion can suffer from Gibbs phenomena (the degree depends to some extent on the diffusivity). Xue (2000) showed that a simple flux limiter could be used in order to ensure a high degree of monotonicity (i.e. greatly diminish over-shoot or Gibbs phenomena) by simply imposing the constraint

$$F^{o-1}_{\varphi\,x} = F^{o-1}_{\varphi\,x} \, {\rm max}\left[ 0,... ..., F^{o-1}_{\varphi\,x} \, \delta_x\varphi\right) \right] \hskip.25in (o \geq 4)$$ (103)

which insures that the diffusion is down-gradient. This rather simple flux limiter gives results consistent which more elaborate schemes at a reduced cost. An example is shown in Fig. 3 for the 6th order diffusion equation of a box function:

$${\frac {\partial Y}{\partial t}} = \kappa^6 {\frac {\partial^6 Y}{\partial x^6}}$$ (104)

As can be seen in Fig. 3, the monotonic flux corrector is quite efficient at removing oscillations/overshoot, at a relatively low cost (and it is relatively simple both conceptually, and code wise).

Figure 3: Solution of the 1-d horizontal diffusion equation (Eq. 105) after a time $t$ using the standard 6th order (red curve) and 6th order using the monotonic limiter from Eq. 104 (blue curve).