Horizontal advection discretization: flux form

ASP contains option for high order horizontal advection schemes which are described herein. The standard second order flux for the horizontal advection of the scalar variable phi defined at a mass point (one-dimensional) is expressed on the C-grid as

$$U_i {\overline\varphi}_i \vert_2 = U_i \left( \varphi_{i}+\varphi_{i-1}\right)/2$$ (134)

where the subscript 2 refers to the order. Note that while simple, horizontal gradients are damped over time and this scheme must be accompanied by additional numerical diffusion (numerical or parameterized).

The main advantage of higher order schemes is the capability to resolve or maintain higher horizontal gradients. The fourth and sixth order advection fluxes can be expressed following Wicker and Skamarock (2002), respectively, as

$$U_i {\overline\varphi}_i \vert_4 = U_i \left[ 7 \left( \varphi_{i} +\varphi_{i-1}\right) - \left( \varphi_{i+1}+\varphi_{i-2}\right) \right]/12$$ (135)

$$U_i {\overline\varphi}_i \vert_6 = U_i \left[ 37\left( \varphi_{i... ...+1}+\varphi_{i-2}\right) + \left( \varphi_{i+2}+\varphi_{i-3}\right) \right]/60$$ (136)

These higher order stencils consist in an economical way to model higher order advection compared to more complex methods (which include interpolating $U$ also). They represent a reasonable compromise between simplicity (and computational speed: it represents a negligible CPU increase compared to the second order flux) and accuracy (results are still comparable to more complex high order stencils).

The above expressions can be used to derive odd-ordered advection schemes. A fifth order scheme can be expressed simply as the sixth order advection scheme less a fourth-order dissipative component Wicker and Skamarock (2002) as

$$\begin{split} U_i {\overline\varphi}_i \vert_5 =& U_i \left[ ... ...left( \varphi_{i+3}-\varphi_{i-2}\right) \right]/60 \end{split}$$ (137)

Note that if 5th order advection is in use, horizontal diffusion (numerical) can be OFF as this scheme is diffusive. One may or may not wish to use (user option) a physically-based or a high order horizontal diffusion/turbulence scheme with this advection scheme (the use of such a diffusion scheme might be of more interest for relatively high spatial resolution runs). Also, note that the diffusive component in Eq. 138 is proportional to the wind speed, so that in low wind speed regions it might be of interest to have some high order or weak low-order horizontal diffusion for numerical (noise suppression) reasons. For example, J. C. Knievel and Hacker (2007) showed that using the aforementioned 5th order scheme together with 6th order horizontal diffusion reduced noisiness in low wind speed regions without impacting the fields much in other regions (for more information on horizontal diffusion, see Section 2.2.2). Indeed, with a reasonable choice for the value of the diffusivity, the 6th order diffusion essentially acts as a noise filter for low wavelength "noise" features (with relatively little damping at higher wavelengths): so that gradients are left fairly well in tact. Therefore currently in ASP, the default is to use 6th order monotonic diffusion (with a low diffusivity) together with the 5th order advection scheme.

Note that one drawback of using higher order schemes is oscillatory behavior (Gibbs phenomena), especially in regions of strong gradients. This is especially problematic for conserved scalars: negative values can result from the above schemes. In ASP, this is corrected using a simple flux limiter (from MesoNH), however, oscillations can still be pronounced at times (as the flux limiter schemes only remove values below some minimum: spiky maximums can still occur). To combat this problem, a more recent class of non-oscillatory schemes has been developed. In ASP, there is an additional high order option which has been implemented for scalars): the 5th order WENO (Weighted Essentially Non-Oscillatory: Liu and Chan (1994)) scheme.

The WENO scheme is an upwind scheme, thus we define the fluxes at the mass cell edges (the stencil depends on the flow direction):

$$U_i {\overline\varphi}_i\vert_5 = U_i {\overline\varphi}_i \left( F^-_{i+1/2}\vert_5 , F^+_{i+1/2}\vert_5 \right)$$ (138)

where $F^-_{i+1/2}$ represents the flux coming from the left (i.e. $U \geq 0$). So, we can also express Eq.139 as

$$U_i {\overline\varphi}_i\vert_5 = U_i {\overline\varphi}_i \left[... ...a_F F^-_{i+1/2}\vert_5 \,+\, \left(1-\delta_F\right) F^+_{i+1/2}\vert_5 \right]$$ (139)

where $\delta_F=1$ if $U \geq 0$, otherwise $\delta_F=0$. The left extrapolated 5th order flux at $F^-_{i+1/2}$ is then defined as a weighted sum of 3 polynomials

$$U^-_{i+1} {\overline\varphi}_{i+1}\vert_5 = U_{i+1} \sum_{l=0}^{2} \omega_l \, v_l$$ (140)

where the polynomials are defined as

$$v_0 = \left( 2 \varphi_{i} + 5 \varphi_{i+1} - \varphi_{i+2} \right)/6$$ (141)

$$v_1 = \left( -\varphi_{i-1} + 5 \varphi_{i} + 2 \varphi_{i+1} \right)/6$$ (142)

$$v_2 = \left( 2 \varphi_{i-2} - 7 \varphi_{i-1} + 11 \varphi_{i} \right)/6$$ (143)

As opposed to the stencils described before the WENO scheme, here the stencil polynomial weights, $\omega$, vary in time and space in such a sway as to minimize oscillations. The WENO weights are defined as

$$\omega_l = {\frac {\alpha_l}{\left(\sum_{l=0}^{2} \alpha_l \right) }}$$ (144)

where

$$\alpha_0 = {\frac{3}{ 10{\left(\epsilon + \beta_0\right)}^2} }$$ (145)

$$\alpha_1 = {\frac{6}{ 10{\left(\epsilon + \beta_1\right)}^2} }$$ (146)

$$\alpha_2 = {\frac{1}{ 10{\left(\epsilon + \beta_2\right)}^2} }$$ (147)

and $\epsilon$ is a small numerical value (to avoid division by zero). Finally, the so-called smoothness indicators, $\beta$, are defined as

$$\beta_0 = {\frac{13}{12}} {\left( \varphi_{i} - 2\varphi_{i+1} - \varphi_... ...{\frac{1}{4}} {\left( 3 \varphi_{i} - 4\varphi_{i+1} + \varphi_{i+2} \right)}^2$$ (148)

$$\beta_1 = {\frac{13}{12}} {\left( \varphi_{i-1} - 2\varphi_{i} - \varphi_{i+1} \right)}^2 + {\frac{1}{4}} {\left( \varphi_{i-1} - \varphi_{i+1} \right)}^2$$ (149)

$$\beta_2 = {\frac{13}{12}} {\left( \varphi_{i-2} - 2 \varphi_{i-1} - \varp... ...+ {\frac{1}{4}} {\left( \varphi_{i-2} - 4\varphi_{i-1} + \varphi_{i} \right)}^2$$ (150)

As the fluxes are upwind in this scheme, the right extrapolated flux, $F^+_{i+1/2}\vert_5 = U^+_{i+1} {\overline\varphi}_{i+1}\vert_5$, is obtained by shifting the stencil by $i+1$ (i.e. to the right) and using symmetry (i.e. by reversing the indices for $\varphi_l$ in the $v_l$ and $\beta_l$ terms.)

A simple example academic test comparing the two 5th order advection schemes presented herein is shown in Fig. 6 and Fig. 7. The 1-dimensional wave equation

$${\frac {\partial Y}{\partial t}} \,+\, u \, {\frac {\partial Y}{\partial x}} = 0$$ (151)

with a constant $u$ value has been integrated (288 time steps) using the (herein named) standard 5th order scheme from Eq. 138. The solution after at time $t$ is shown in red, while the solution using the aforementioned simple flux limiter is shown in green. The oscillatory behavior is well illustrated in this problem (owing to the sharp gradients).

Figure 6: Solution of the 1-d wave equation (Eq. 152) using the standard 5th order advection scheme from Eq. 138. The solution after at time $t$ is shown in red, while the solution using the simple flux limiter is shown in green.

Using the 5th order WENO scheme (Eq. 141), the solution after at time $t$ is shown in red, while the solution using the standard 5th order scheme with the aforementioned simple flux limiter is shown, once again, in green. The non-oscillatory behavior is clearly seen in the WENO solution, although the gradients are slightly less. Also note (not shown) that as $t$ increases, as should be expected, the WENO scheme preserves the wave form much better than the standard scheme. Note that this basic example is rather extreme (box function), as real meteorological phenomena are more likely to have a more smooth character, however, it is a classic numerical problem that highlights the advection scheme differences.

Figure 7: As in Fig. 6, except using the 5th order WENO scheme (Eq. 141) (in red).

In summary, the 5th order WENO scheme seems to be advantageous for conserved scalars compared to the standard 5th order scheme (especially for cloud water, and tracers which might have zero values over considerable regions of the domain and be characterized by potentially very sharp gradients). Also, the flux limiter is technically not needed for the WENO scheme (although it can be applied for numerical, round-off reasons). However, as the WENO scheme requires more computations, continuous fields like the wind components and even temperature could likely still be advected using the standard scheme. In addition, the diffusive behavior of the standard scheme (with diffusion proportional to the wind speed) seems to be advantageous for the main variables of the dynamic core ($T$, $u$, $v$, and also $w$ and $\Phi$ for the non-hydrostatic module) since this diffusion seems to be very efficient at damping any numerical noise which can develop (while still maintaining relatively strong gradients and negating the need for additional or parameterized numerical diffusion). So, currently the default model configuration uses the WENO 5th order scheme for all scalars, and the standard 5th order horizontal advection scheme for the variables controling the dynamics.

As a final note, tests have indicated that to save CPUs, the standard 5th order scheme can be used for scalars for the first 2 RK3 steps, and the 5th order WENO scheme can be used just for the last RK3 step. This configuration has been found to give results quite close to the configuration with the WENO scheme applied at all 3 RK3 steps. Currently, this is the default configuration in ASP.

Finally, the odd-ordered advection schemes can be used owing to the time discretization method utilized by the model (3rd or 4th order Runge-Kutta schemes: see details on time integration in Section 2.10).