Estimation of energy dissipation caused by odd order difference schemes for an unstable planetary boundary layer

Numerical error by an odd order difference approximation in the advection term dissipates energy in the resolved motion. We propose a method to estimate the energy dissipation rate derived from the numerical error to clarify its effect on the energy budget. The energy dissipation rate in Large‐Eddy Simulations for an ideal unstable boundary layer is analysed. The proposed method reveals that the imbalance of the energy budget can be ascribed to the non‐conservative property with respect to kinetic energy of the discretised advection terms, while the net energy dissipation in the resolved field is almost insensitive to the origins of the energy dissipation. However, the magnitude of the momentum and temperature fluxes is sensitive to the difference scheme; the numerical error in the advection term reduces these fluxes in the resolved part. A subgrid scheme should be designed so as to compensate decrease of the resolved flux when an odd order difference is applied to the advection term.


| INTRODUCTION
Large-Eddy simulation (LES) has been widely performed for numerical simulations in the planetary boundary layer due to the advance in the computational power. The momentum and heat fluxes in the unresolved subgrid scale (SGS) are parameterized in the LES, while the resolved grid scale (GS) motions are explicitly solved. Furthermore, the LES provides the same dissipation rates as that experienced by the flow at the GS cutoff scale.
On the other hand, the effect of a numerical error derived from the discretization of the advection term cannot be neglected when finite differences are used in a numerical simulation with a finite number of grid points. The numerical error in an odd order difference scheme possesses a diffusive property, which results in decay in small scales of motion. It causes energy dissipation in the GS motions. The energy is not generally conserved also in even order difference without any diffusive error unless the energy is treated as a prognostic variable in a numerical model, although Morinishi et al. (1998) found a specific form of even order centre difference schemes which simultaneously conserve the momentum and energy. Ghosal (1996) and Chow and Moin (2003) analysed numerical errors in LES in terms of energy spectra. Margolin et al. (1999) estimated the energy conservation error from the imbalance in the energy budget equation. However, an LES parameterization scheme is generally formulated without accounting for numerical error effects.
Since the effect of the energy dissipation due to the numerical error is qualitatively similar to that of an SGS parameterization, it could be used as an alternative to a turbulence scheme without any explicit formulation for the SGS motion. This approach is called implicit LES (ILES). The ability of the ILES is reported in simulations of thermal convection (Margolin et al., 1999;Piotrowski et al., 2009) and stratocumulus clouds (Pressel et al., 2017). Nevertheless, quantitative evaluation of the distribution of energy dissipation from the numerical error is difficult, while the averaged value over an entire domain can be accurately calculated as we discuss later. Additionally, whether the momentum and heat fluxes obtained from ILES are identical to those of an SGS scheme is uncertain.
In the present study, we propose a method for estimating the energy dissipation rate associated with the numerical error in the advection term. Numerical experiments for an ideal unstable boundary layer are performed to evaluate the energy dissipation rate including the effect of the numerical error for various difference schemes. We focus on the energy conservation error in an odd order difference scheme for the advection term. The odd order difference scheme is widely used in atmospheric models like the weather research and forecasting (WRF) model (Wicker and Skamarock, 2002). Furthermore, we analyse the dependence of the momentum and temperature fluxes at the surface and these vertical transfers on the finite difference schemes and discuss the validity of the ILES concept in terms of such transfers.
2 | ANALYSIS PROCEDURE 2.1 | Estimation of the energy dissipation rate due to the numerical error In the present analysis, we consider the full compressive non-hydrostatic equations. It is assumed that the density ρ and the momentum ρu i (u i : the velocity in the i-th direction) are treated as the prognostic variables, and these advection terms are evaluated as a flux form in a numerical model as shown in Equations A1-A4. The governing equation for the GS part of the kinetic energy per unit volume is written as P s , P b and Ξ in Equation 1 represent the shear production term for the SGS turbulent kinetic energy (TKE), the buoyancy production term and the energy transfer written in the flux form, respectively. The overbar denotes the volume average over a grid cell. These definitions are given in Section B, Supporting Information, Data S1. The term ϵ adv , which is mathematically zero, are added for the following reason.
The first and second terms of ϵ adv defined in Equation 2 correspond to the tendency of the kinetic energy due to the advection terms evaluated in a numerical model based on Equations A1-A4. Note that the second term is zero if the fluid is assumed to be incompressible. The third term is the flux form of the first and second terms. The term ϵ adv in a discretised form is not zero due to the discretised error, although it should be zero because the sum of the first and second terms is equivalent to the third term in the continuous form Equation 2. Nevertheless, the volume average of the third term is also zero in the discretised form without the energy flux from the boundary, and the volume average of ϵ adv is not affected by the numerical error of the third term. If the volume average of ϵ adv is not zero, this error can be attributed to the first and second terms. ϵ adv can be interpreted as the net dissipation rate of the kinetic energy excluding the energy flux which does not change the total energy. In the present study, we define ϵ adv as the dissipation rate derived from the numerical error of the advection term under the above assumption. However, ϵ adv should be regarded as an estimated value for the dissipation rate, because its profile is distorted by the numerical error of the third term which cannot be excluded. We discuss the horizontal average of the energy dissipation in the following section. The horizontal averaged error from the third term of ϵ adv appears only in the vertical difference, because the terms of the horizontal difference written in a flux form vanish by the integration over a horizontal plane. Note that ϵ adv can be divided into the contributions from the horizontal and vertical derivatives ϵ adv = ϵ adv(hor.) + ϵ adv(ver.) because the right-hand side in Equation 2 is mathematically zero for each j component.

| Experimental setup
We perform the numerical experiments for an unstable boundary layer to evaluate the energy dissipation due to the error of the discretised advection term in an LES. The numerical model employed in the experiments is ASUCA, which was developed and is operated as a regional numerical weather prediction (NWP) system by Japan Meteorological Agency (JMA) (Ishida et al., 2010;Aranami et al., 2015). The dynamical core of ASUCA is briefly described in Section A in Data S1. The advection term in ASUCA is evaluated by the third upwind difference with the flux limiter by Koren (1993) in operational use.
In the present study, the Kawamura-Kuwahara scheme (Kawamura and Kuwahara, 1984) (hereafter referred to as the KK scheme) with the third and fifth orders is employed to calculate the advection term because the scheme allows us to systematically control the coefficient in the leading term of the error in a differential approximation for the spatial derivative. When the spatial derivative at the i-th stencil is approximated with n-th order accurate as the KK scheme approximates ϕ i + 1/2 as represents the fourth-order centre difference. The non-dimensional parameter α is introduced in this scheme. The KK scheme is identical to the third-order upwind difference using three stencils when α = 1 and the fourth-order centre difference when α = 0. The leading term of the error While the original KK scheme is the third-order accurate, the approximation with the fifth-order can be constituted in a similar manner: The leading term of the error is − αΔ 5 x =60 À Á ∂ 6 ϕ=∂x 6 À Á i . Ideal LES experiments for an unstable boundary layer, whose configurations follow Kitamura and Ito (2016), are examined. The horizontal mean of the initial temperature is given as Θ(z) = 288.15 + 4 × 10 −3 z (z: the height) to impose stable stratification. The geostrophic wind speed U g is set to 4m/s and the initial velocity is identical to the geostrophic velocity in the entire domain. The Coriolis parameter f is 1.031 × 10 −4 /s. The surface temperature is set as Θ s = 291.15 + 20sin(πt/12) K (t: the elapsed time in hours). The surface flux is diagnosed based on the stability function by Beljaars and Holtslag (1991). The roughness lengths for the momentum and heat are assumed to be 0.01 and 0.001 m, respectively. The SGS model proposed by Brown et al. (1994) is employed. The lateral boundary condition is doubly cyclic so that there are no net mass and momentum fluxes from the lateral boundary. The uniform grid size is 20 m in all directions. The numbers of grid points are 480 and 160 in the horizontal and vertical directions, respectively. A sponge layer is placed in the upper 20% of the domain to avoid artificial reflection of gravity waves at the top of the boundary. Time integration is performed for 5 hr with the third-order Runge-Kutta scheme and the result for the last 1 hr obtained with every 1 min is analysed. The time interval for the numerical integration is 1.5 s.; the Courant number is 0.3 for the geostrophic wind.

| RESULTS
The upper panels of Figure 1 show the horizontal average of the budget evaluated as the discretised form of the GS kinetic equation Equation 1 and the estimated numerical error of the advection term ϵ adv . The obtained imbalance in the energy budget P b − P s + Ξ is consistent with Margolin et al. (1999) and close to the error derived from the advection term ϵ adv for all the experiments. This result suggests that the imbalance in the energy budget can be ascribed to the non-conservative property with respect to kinetic energy of the discretised advection terms and ϵ adv successfully extracts such the aspect in the present experiments. The effects of ϵ adv are not negligible in comparison with the energy transfer to the SGS TKE due to the shear production P s . In particular, the numerical error is much larger than the SGS term in the ordinary upwind difference schemes (α = 1), suggesting the importance of considering the effect of the numerical error as well as the SGS parametrization when using an odd order difference scheme. It is confirmed that the energy conservation error in the fourth (sixth) order is smaller than that in the third (fifth) order by two order of magnitude except in the vicinity of the ground (not shown). The non-zero value of ϵ adv almost comes from the term constituting the odd order difference in Equation 4 or 6.
Each term of ϵ adv is plotted in the lower panels of Figure 1. ϵ adv decreases with decreasing α for the same order of difference and ϵ adv with the fifth-order is smaller than that with the third-order for the same α. This result is consistent with the magnitude of the leading error term. The second term in Equation 2 is much smaller than the other terms; the effect of the compressibility is negligible in the present experiments.
The dependence of the energy dissipation rate on α is plotted in Figure 2. The averaged value from z = 100 m to the height of the boundary layer defined as the level of the minimum total temperature flux is displayed. Note that the estimate under z = 100 m is excluded in the analysis because the accuracy of the estimation of the numerical error appears to be uncertain in this region. The vertical part of the error ϵ adv(ver.) , in which the numerical error in the flux term (the third term in ϵ adv ) cannot be excluded, is sufficiently smaller than the horizontal part for all experiments, suggesting that ϵ adv can approximate the actual dissipation rate in the mixed layer. This figure clearly shows that the net energy dissipation rate is insensitive to the parameter α and the order of the spatial difference. We examine an additional experiment in the pure ILES without any SGS model. In this experiment, the energy is completely dissipated by the numerical error. The energy dissipation rate in the pure ILES is similar to the result including the SGS scheme as shown by the dashed line in Figure 2. The energy dissipation of the numerical error can serve as an alternative to the energy transfer to the SGS component represented by an SGS scheme, and the ILES concept can be justified in terms of the budget of the resolved energy. Although the net energy dissipation is almost independent of a spatial difference scheme, the obtained GS motions are not necessarily similar to each other. The spectrum of the kinetic energy at z = 500 m is shown in Figure 3. The energy dissipation range in higher wavenumbers is expanded, and its slope is steeper with increasing α. The fifth-order scheme produces a more apparent inertial subrange with −5/3 slope than the third-order scheme with same α. This result is consistent with Mittal and Moin (1997). It implies that the dissipation due to the numerical error of the advection term works in wider wavenumbers than that due to the SGS parameterization. The wider dissipation range in the wavenumber space causes degradation of the effective resolution in a numerical simulation.
Next, we analyse the dependence of the momentum and temperature fluxes on a spatial difference scheme. Figure 4 displays the momentum and temperature fluxes at the surface. The absolute value of both surface fluxes decreases with increasing α and the flux in the third-order scheme is smaller than that in the fifth-order scheme; the surface flux decreases as ϵ adv becomes larger. The spread of the surface flux in the experiments is approximately 15% for the momentum and 4% for the temperature. In contrast to the insensitive net energy dissipation to the advection scheme, the momentum flux is especially affected by the selection of the advection scheme. Matheou (2016) reported that the effect of the order of accuracy could significantly affect the surface flux in an LES for a stable boundary layer. The vertical profiles of the temperature and its vertical flux simulated with the fifth-order difference scheme are shown in Figure 5. The potential temperature in the mixed layer becomes lower for a larger α because of decreased the heat supply from the surface. The GS part of the temperature flux is more sensitive to the numerical error than the SGS part, especially in the surface layer and the top of the boundary layer. The GS flux in the surface layer is smaller for the α = 1 case. The degradation of the effective resolution observed in Figure 3 may be attributed to a decreased GS flux.
On the other hand, decreasing the production of the SGS TKE P s by increasing ϵ adv results in a decreased SGS TKE and eddy viscosity coefficient K m ; thus contributing to a decreased SGS flux. Therefore, the GS and SGS fluxes are not mutually complementary, although P s can compensate for ϵ adv with respect to the GS energy budget. The result suggests that the ILES concept is not valid for evaluating the momentum and temperature fluxes and the SGS parameterization must be carefully designed according to a difference scheme applied to the advection term; an SGS scheme should be designed so as to compensate decrease of the GS flux. Matheou (2016) indicated that the surface fluxes depend on model constants for an SGS scheme in the experiments of a stable boundary layer; the lower Smagorinsky constant leads to over-prediction of the surface heat flux when an even order difference scheme is employed. Notably, the temperature flux does not converge to zero above the boundary layer for the α = 0.1 experiment due to the insufficient energy dissipation. The SGS scheme by Brown et al. (1994) reduces the eddy viscosity coefficient as static stability is stronger because of negative buoyancy production in the SGS TKE budget. It might be speculated that the GS energy dissipation becomes insufficient and artificial disturbance emerges above the boundary layer without the dissipation due to the numerical error.

| CONCLUDING REMARKS
In the present study, we propose a method to estimate the energy dissipation rate derived from the numerical error appearing in the advection term approximated by an odd order difference. Ideal experiments for an unstable boundary layer are performed with the KK scheme to evaluate the advection term, and the horizontal mean of the energy dissipation rate is analysed. The proposed method seems to adequately estimate the effect of the numerical error on the energy dissipation at least for an ideal convective boundary layer. The estimated energy dissipation due to the numerical error ϵ adv increases as α increases and compensates for this increase by the SGS scheme P s ; the net energy dissipation in the GS motion is insensitive to the difference scheme in the boundary layer. The ILES concept, in which the dissipative property in the numerical error is regarded as an SGS parameterization, can be justified in this sense.
On the other hand, the spreads of the momentum and temperature fluxes are about 15 and 4%. In particular, the effect of the numerical error in an odd order scheme is not negligible for the momentum flux. The GS part of them is smaller when the energy dissipation from the numerical error is larger, possibly because the numerical error dissipates the energy in a wider wavenumber range than the SGS scheme. The SGS flux does not function to compensate for a shortage of the GS flux; this problem is inevitable if an odd order scheme is applied. At least the effect of a difference scheme used in the advection term should be taken into account in the design of an SGS scheme to compensate decrease of the GS flux.
Although selecting a small α in the KK scheme can alleviate the effect of the energy dissipation by the numerical error in an odd order accuracy difference, it also facilitates the artificial disturbance frequently observed in an even order difference scheme. An additional explicit (hyper-)viscosity term is often introduced when an even order difference scheme is employed (e.g., Nishizawa et al., 2015). However, introducing the explicit hyper-viscosity cannot solve the consistency in the flux because the viscosity term controls the energy dissipation in the GS motion and the appropriate value of the viscosity coefficient must be determined empirically. The length scale in a LES would be an effective resolution rather than a grid size itself as Lilly (1967) discussed. While this idea suggests that a larger length scale should be employed when an odd order difference scheme is used, the effective resolution corresponding to a spatial difference scheme would be difficult to estimate. A further study is necessary to explore how to appropriately determine the length scale. Wicker, L.J. and Skamarock, W.C. (2002)

SUPPORTING INFORMATION
Additional supporting information may be found online in the Supporting Information section at the end of this article.