[MGI40-P01] A study on numerical accuracy of atmospheric dynamical cores necessary for large eddy simulations
Keywords:atmospheric planetary boundary layer turbulence, large eddy simulation, numerical accuracy of discretization for dynamical core
1. Introduction
As computational resources of supercomputers increase, spatial scale which global atmospheric models can resolve has been recently approaching that in large eddy simulations (LES). Then, the authors consider that one of issue is the accuracy of numerical schemes applied for the dynamical core of nonhydrostatic atmospheric models. Because it is possible that the numerical dissipation associated with advection schemes by low-order finite volume method dominates eddy mixing parametrized by the Smagorinsky-Lilly type turbulent model. We derive a criterion of numerical accuracy necessary for LES of atmospheric planetary boundary layer turbulence, based on the theory of 3-dimensional isotropic homogeneous turbulent flow, and check whether the criterion is reasonable in actual simulations.
2. Criterion for numerical accuracy necessary to LES
For the spatial scale, lΔx (where l=1,2,.., and Δx is grid spacing), the ratio of e-folding time of decay due to eddy viscosity term to that due to numerical diffusion term with 2n differential operator is derived:
R= (m/π)4/3 * (2πCs)2 /(|U| γadv) * (l/π)2(n-1) * (Δx)1/3
where the magnitude of strain velocity tensor in eddy viscosity is related to Δx, following the theory of 3-dimensional isotropic homogeneous turbulent flow (note that we determine the value of η based on the results from numerical experiments), Cs is the Smagorinsky constant, m is the ration of spatial filter length to Δx and γadv is non-dimensional strength of numerical diffusion related to the background advective velocity, |U|. For upwind schemes, the numerical diffusion is implicitly determined proportional to wind velocity, while for central schemes it is explicitly added. In order for the numerical diffusion not to dominant the eddy viscosity, It is important that R>>1 in effective resolution. However, the R tends to reduce as spatial resolution increases. Next, we investigate the influence of numerical diffusion on calculations of typical atmospheric LES, and consider the range of R realized.
3. Investigation by numerical experiment
[Experimental setup]
We perform idealized numerical experiments same as ones in Nishizawa et al. (2015), in which we use both of central and upwind advection schemes with several spatial accuracy. A regional nonhydrostatic model, SCALE-RM, is used in this study. The size of computational domain is 9.6 x 9.6 x 3 km3. As the initial condition, potential temperature disturbance is added to stably stratified atmosphere. The initial background horizontal wind is set to 5 m/s, and heat flux with 200 W/m2 is continuously given at surface. As advection schemes, 3rd and 5th order upwind schemes (UD3, UD5), and 2nd, 4th and 6th order central schemes (CD2, CD4 and CD6) are used. For central schemes, numerical diffusion is explicitly added. The e-folding time is about 1 minute for 2-grid scale, and the number of differential operater (ND) is basically set to 2, 4 and 6 so that leading dispersive error terms are eliminated for CD2, CD4 and CD6, respectively. We also conduct a experiment with CD4ND8 as in Nishizawa et al. (2015). For all cases, the grid size is 10m, and the temporal integrations are performed over 4 hours.
[Result]
Figure (a) shows that, in irrespective of advection schemes, the energy spectra in inertial subrange are lower than that following -5/3 power law at wavelength shorter than 8-10 times grid spacing (referred to as “effective resolution” here). Furthermore, due to the implicit diffusion proportional to background wind speed, the fields for upwind schemes are more diffusive than that for central schemes. Over the wavelength shorter than effective resolution, the energy spectra for 3rd and 5th order upwind schemes are lower than 2nd and 4th order central schemes with weak numerical filter, respectively. Our examination suggests, if upwind scheme is used, it is desirable that the accuracy of spatial discretization for advection terms is 5th order or higher. Figure (b) shows the dependence of R on Δx and l, which is consistent to the result from numerical experiment mentioned above. It indicates that R we derived here provides a guide for numerical accuracy of advection schemes necessary for LES.
4. Future work
We will investigate the influence of temporal discretization, and will evaluate the effect of “perfectly” high-order discretization where all terms are spatially and temporally discretized with accuracy precisely higher than 2nd order. To do so, we will also explore suitable numerical schemes for dynamical cores.
As computational resources of supercomputers increase, spatial scale which global atmospheric models can resolve has been recently approaching that in large eddy simulations (LES). Then, the authors consider that one of issue is the accuracy of numerical schemes applied for the dynamical core of nonhydrostatic atmospheric models. Because it is possible that the numerical dissipation associated with advection schemes by low-order finite volume method dominates eddy mixing parametrized by the Smagorinsky-Lilly type turbulent model. We derive a criterion of numerical accuracy necessary for LES of atmospheric planetary boundary layer turbulence, based on the theory of 3-dimensional isotropic homogeneous turbulent flow, and check whether the criterion is reasonable in actual simulations.
2. Criterion for numerical accuracy necessary to LES
For the spatial scale, lΔx (where l=1,2,.., and Δx is grid spacing), the ratio of e-folding time of decay due to eddy viscosity term to that due to numerical diffusion term with 2n differential operator is derived:
R= (m/π)4/3 * (2πCs)2 /(|U| γadv) * (l/π)2(n-1) * (Δx)1/3
where the magnitude of strain velocity tensor in eddy viscosity is related to Δx, following the theory of 3-dimensional isotropic homogeneous turbulent flow (note that we determine the value of η based on the results from numerical experiments), Cs is the Smagorinsky constant, m is the ration of spatial filter length to Δx and γadv is non-dimensional strength of numerical diffusion related to the background advective velocity, |U|. For upwind schemes, the numerical diffusion is implicitly determined proportional to wind velocity, while for central schemes it is explicitly added. In order for the numerical diffusion not to dominant the eddy viscosity, It is important that R>>1 in effective resolution. However, the R tends to reduce as spatial resolution increases. Next, we investigate the influence of numerical diffusion on calculations of typical atmospheric LES, and consider the range of R realized.
3. Investigation by numerical experiment
[Experimental setup]
We perform idealized numerical experiments same as ones in Nishizawa et al. (2015), in which we use both of central and upwind advection schemes with several spatial accuracy. A regional nonhydrostatic model, SCALE-RM, is used in this study. The size of computational domain is 9.6 x 9.6 x 3 km3. As the initial condition, potential temperature disturbance is added to stably stratified atmosphere. The initial background horizontal wind is set to 5 m/s, and heat flux with 200 W/m2 is continuously given at surface. As advection schemes, 3rd and 5th order upwind schemes (UD3, UD5), and 2nd, 4th and 6th order central schemes (CD2, CD4 and CD6) are used. For central schemes, numerical diffusion is explicitly added. The e-folding time is about 1 minute for 2-grid scale, and the number of differential operater (ND) is basically set to 2, 4 and 6 so that leading dispersive error terms are eliminated for CD2, CD4 and CD6, respectively. We also conduct a experiment with CD4ND8 as in Nishizawa et al. (2015). For all cases, the grid size is 10m, and the temporal integrations are performed over 4 hours.
[Result]
Figure (a) shows that, in irrespective of advection schemes, the energy spectra in inertial subrange are lower than that following -5/3 power law at wavelength shorter than 8-10 times grid spacing (referred to as “effective resolution” here). Furthermore, due to the implicit diffusion proportional to background wind speed, the fields for upwind schemes are more diffusive than that for central schemes. Over the wavelength shorter than effective resolution, the energy spectra for 3rd and 5th order upwind schemes are lower than 2nd and 4th order central schemes with weak numerical filter, respectively. Our examination suggests, if upwind scheme is used, it is desirable that the accuracy of spatial discretization for advection terms is 5th order or higher. Figure (b) shows the dependence of R on Δx and l, which is consistent to the result from numerical experiment mentioned above. It indicates that R we derived here provides a guide for numerical accuracy of advection schemes necessary for LES.
4. Future work
We will investigate the influence of temporal discretization, and will evaluate the effect of “perfectly” high-order discretization where all terms are spatially and temporally discretized with accuracy precisely higher than 2nd order. To do so, we will also explore suitable numerical schemes for dynamical cores.