Hamiltonian tools applied to non-hydrostatic modeling

1. Hamiltonian tools applied to non-hydrostatic modeling Almut Gassmann Max Planck Institute for Meteorology Hamburg, Germany

2. Motivation... Max Planck Institute for Meteorology and German Weather Service (DWD) are developing a new global model system in a joint project: ICON Members and collaborators present at IPAM: Marco Giorgetta, Peter Korn, Luis Kornblueh, Leonidas Linardakis, Stephan Lorenz, Almut Gassmann, Werner Bauer, Florian Rauser, Hui Wan, Peter Dueben, Tobias Hundertmark, Luca Bonaventura My part: non-hydrostatic atmospheric model Contents... Hamiltonian form for the continuous moist turbulent equations Spatial discretisation of Poisson brackets Temporal discretisation Contents • Hamiltonian form • Spatial discretisation • Temporal discretisation

3. Governing equations for NWP and climate simulations Is there really a need to reconsider that? No! But: As soon as moisture and turbulence averaging come into play, things become ugly: → Approximations to thermodynamics may distort local mass and/or energy consistency. Unfortunately, we do not know the longer term impact of such small, but systematic errors. Contents • Hamiltonian form • Spatial discretisation • Temporal discretisation

4. Starting from general equation set... Momentum equation Conservation of total (moist) mass First law of thermodynamics Conservation of tracer mass q: specific quantities Approximations that do not change mass or energy balance: - neglect the molecular heat flux against the turbulent one: W → R - neglect the molecular dissipation against the turbulent one Make sure that the diffusion fluxes of the constituents and the conversion terms sum up to 0. Contents • Hamiltonian form • Spatial discretisation • Temporal discretisation

5. The energy budget must be closed ... shear production buoyancy production dissipation rudimentary mean turbulent kinetic energy equation This suggests: adding in the heat equation Contents • Hamiltonian form • Spatial discretisation • Temporal discretisation

6. Specifics for a moist atmosphere... The internal energy u is not a suitable variable. We have to unveil the phase changes of water or chemical reactions, and thus consider rather the enthalpy h instead of the internal energy u. The ideal equation of state is assumed to be valid also for averaged quantities virtual increment Finally, a prognostic temperature equation is obtained. Contents • Hamiltonian form • Spatial discretisation • Temporal discretisation

7. Towards flux form equations... viewpoints of fluid dynamics • Particle view (Lagrange) • Field view (Euler) • Flux form for scalars • specific moisture quantities • some form of entropy variable • Because the actual entropy s including all moisture quantities is impractical to handle, we decide for a compromise: the virtual potential temperature. • We write the wind advection in Lamb form to unveil the vorticity (reason: particle relabelling symmetry). The views are equally valid and suitable for building a numerical model. In our ICON project, we decide for the Eulerian standpoint. Contents • Hamiltonian form • Spatial discretisation • Temporal discretisation

8. Equation set that unveils the entropy production... „entropy production“ virtual potential temparture Next step: Poisson bracket form for the non-dissipative adiabatic limit case... Contents • Hamiltonian form • Spatial discretisation • Temporal discretisation

9. Hamiltonian dynamics ... A suitable Hamiltonian functional at least covers the adiabatic part of the dynamics: The Hamiltonian is a function of the density, a suitable thermodynamic variable and the velocity. With the choice of the virtual potential temperature density we obtain ''dynamic'' and ''thermodynamic'' functional derivatives independently. Hamiltonian dynamics (at least for the adiabatic part) Contents • Hamiltonian form • Spatial discretisation • Temporal discretisation

10. Poisson bracket... scalar triple product is antisymmetric: A.(BxC) = -B.(AxC) = -C.(BxA) Antisymmetry: swapping F and H only alters the sign, F=H gives a zero bracket result Note: only the divergence operator appears, not the gradient operator, duality of the div and grad operators is automatically given Background: integration by parts rule Contents • Hamiltonian form • Spatial discretisation • Temporal discretisation

11. Summary so far... • Exact Hamiltonian form with Poission brackets seems to be only practicable to write down for idealized (dry, non-turbuent) flows. • This is no contradiction to the conservation of the total energy, because friction contributes correctly also to the internal enery by dissipative heating, and phase changes do not change the total energy. • The structure of the Poisson bracket guarantees for correct energy conversions and thus energy conservation comes as a by product. • Mass conservation is automatically given. • The virtual potential temperature enters the equations as a passive tracer – as expected. In the dry adiabatic non-dissipative limit case, the entropy is conserved. • Prognostic variables might be chosen freely. „Nice“ prognostic variables are the density and virtual potential temperture density (also the Exner pressure). Next step: Discretize brackets instead of single terms in the equations... Contents • Hamiltonian form • Spatial discretisation • Temporal discretisation

12. Numerics with brackets... R. Salmon (2005, etc.): „...From the standpoint of differential equations, conservation laws arise from manipulations that typically include the product rule for derivatives. Unfortunately, the product rule does not generally carry over to discrete systems; try as we might, we will never get digital computers to respect it. However, in the strategy adopted here, conservation laws are converted to antisymmetry properties that transfer easily to the discrete case; digital computers understand antisymmetry as well!“ Contents • Hamiltonian form • Spatial discretisation • Temporal discretisation

13. C-grid discretisations... Poisson brackets convert easily the Arakawa C grid (also in the vertical ➞ Lorenz grid): Requirements: - divergence via Gauss theorem - Laplacian-consistent inner product Note: θv is not touched by the bracket philosophy, it is only required 'somehow' at the interface position. Contents • Hamiltonian form • Spatial discretisation • Temporal discretisation

14. prognostic equations Role of the potential temperature as a 'tracer'... • Higher order advection scheme is interpreted as to give an interface value for theta. • Well balancing approach (for terrain following-coordinates) might be interpreted as to give a specialnearly hydrostatic state via the estimation of theta at the vertical interface. • A combination of both requirements is also possible -> next slide. Contents • Hamiltonian form • Spatial discretisation • Temporal discretisation

15. Well balancing... • Atmospheric motions are nearly hydrostatic. The local truncation error of the pressure gradient term might violate the nearly hydrostatic state. Workaround: locally well balanced reconstruction with the help of a local hydrostatic background state (Botta et al., 2002) contrib. to covariant horizontal equation: interface value in flux divergence term: Contents • Hamiltonian form • Spatial discretisation • Temporal discretisation

16. Vorticity flux term... absolute potential vorticity gives the mass flux • Same type of vector reconstruction. • In case of an (irregular) hexagonal grid, further consistency requirements are required, which determine the stencil and method for the vector reconstruction. • The PV takes the same role as theta in the previous considerations: it is the tracer quantity in the vorticity equation and might be subject to further conditions (anticipated vorticity flux method etc.). Contents • Hamiltonian form • Spatial discretisation • Temporal discretisation

17. The nonlinear advection of momentum is split into two parts. vorticity flux term gradient of the kinetic energy belongs to the 'divergent' part of the flow Time integration scheme... State of the art nonhydrostatic models • fully implicit (expensive and global) • vertically implict • horizontally explicit (forward-backward wave solver in combination with a Runge-Kutta type scheme for advection: split-explicit) How does the time integration scheme look like, if we have discretized Poisson brackets? Contents • Hamiltonian form • Spatial discretisation • Temporal discretisation

18. Energy conserving time integration scheme... Shallow water example Energy budget equation Contents • Hamiltonian form • Spatial discretisation • Temporal discretisation

19. Forward-backward time integration scheme for waves. But kinetic energy term is implicit! Goal: Explicit time stepping... Contents • Hamiltonian form • Spatial discretisation • Temporal discretisation

20. Predictor step philosophy... • similar to RK2 procedure • proposed prediction step Linear implicit method behaves similarly to RK2 split explicit scheme and is unstable. Linear implicit method for the predictor step behaves similarly to the proposed explicit scheme and is stable. Linear stability analysis reveals an unstable behaviour, which is not found in numerical experiments, presumably because the whole scheme is nonlinear. The predictor step plays the role of divergence damping in traditional split-explicit methods. Contents • Hamiltonian form • Spatial discretisation • Temporal discretisation

21. Pressure gradient term... pressure gradient term 'unfortunately implicit' Trick to make it explicit relabeling of time levels gives explicit scheme Remark: This approach is known with empirical weights as acoustic mode filtering (Klemp et al, 2007). We can shed different light on this procedure. Our new weigths are physically based. Contents • Hamiltonian form • Spatial discretisation • Temporal discretisation

22. Summary on the temporal discretisation... • Product rule for derivatives in the context of time integration. • Strict splitting between 'wave solver' and 'advection' becomes questionable. • New light is shed on split-explicit schemes: Alternative explanations for • divergence damping • acoustic mode filtering • The vorticity flux term is still an outsider here. Because it should be energetically neutral, the mass flux therein must be consistent with the continuity equation – also in the time level choice. The time level of the PV itself is not constrained.

23. Baroclinic wave test case Non-hydrostatic ICON model on the hexagonal grid (dx = 240km) including some of the numerical issues disussed in the talk. The run is without additional diffusion.