Researchers are constantly working to further improve the accuracy, efficiency, and scalability of the dynamical core, the physics, and their coupling. The rapid development of computing systems towards massive use of graphics processing units and extreme parallelism requires adaptation of algorithms to further increase model efficiency via strong or weak parallel scaling. Recent years have also seen a rapid increase in hybrid approaches that combine physics-based modelling with data driven techniques, adopting techniques from scientific machine learning for Earth system modelling.
This session invites presentations on the development, testing, and application of novel numerical techniques for Earth system models in a broad sense. The scope includes modifications to the governing equations, horizontal and vertical discretizations, structure preserving methods, time stepping schemes (including parallel in time schemes), advection schemes, adaptive multi-scale models, physics-dynamics coupling, regional and global models, classical and stochastic physical parameterizations, as well as hybrid schemes combining numerical methods and machine learning.
Orals: Tue, 5 May, 08:30–10:15 | Room 1.61/62
The pursuit of global kilometer-scale modeling in atmosphere science presents a dual challenge: comprehensive models become computationally prohibitive and analytically intractable, while idealized models lack essential physical forcing, creating a critical gap in atmospheric research. To bridge this divide, this work develops and validates the innovative ACC-LMARS-SW model, a novel theoretical framework designed for efficient, high-resolution global atmospheric dynamics studies. The model features three synergistic innovations. First, its physics kernel introduces real-world topography and a real-world-equivalent season-aware radiative forcing scheme into the shallow-water equations, significantly enhancing physical realism while retaining conceptual clarity. Second, it employs a variable-resolution stretched cubed-sphere grid that concentrates computational resources on regions of interest , achieving sub-kilometer local resolution without prohibitive global cost. Third, the dynamical core is fully redesigned for massive parallelism using the OpenACC standard, leveraging the efficiency of the Low-Mach Approximate Riemann Solver (LMARS) to harness GPU acceleration.
A suite of numerical experiments demonstrates the model's capabilities. 1) Mesh efficacy: Stretched-grid simulations show reduce error in target areas compared to uniform-resolution runs in classic tests, better resolving nonlinear eddy interactions. 2) Physical fidelity: A global simulation at ~1.5 km average resolution (with ~500 m over the South China Sea) forced by real topography and idealized radiation spontaneously generates a horizontal kinetic energy spectrum featuring distinct k-3 and k-5/3 power-law segments , which is a hallmark of real atmospheric scale interactions. 3) Computational performance: GPU implementation achieves up to a 70x speedup versus estimated serial CPU execution, making global kilometer-scale theoretical experiments feasible on a single workstation.
In conclusion, ACC-LMARS-SW successfully integrates algorithmic design, mesh technology, and high-performance computing to create an efficient and physically insightful "numerical laboratory." This approach provides the high resolution and computational efficiency needed to study fundamental multi-scale dynamics, such as cross-scale energy transfer, dynamical responses to orography, and mesoscale eddy dynamics.
How to cite: Li, Z. and Chen, X.: Accelerated LMARS Shallow Water Model with a Stretched Cubed Sphere Grid at Sub-km Scales, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-4665, https://doi.org/10.5194/egusphere-egu26-4665, 2026.
Radiative transfer plays a key role in the atmosphere’s thermodynamics as photons from the Sun interact with the atmosphere and with other photons being emitted by the earth itself. In climate modeling and numerical weather prediction, it is very common to simplify radiative transfer into a one-dimensional, purely vertical, phenomenon for the sake of saving computational resources. This is because solving the radiative transfer equation requires building and solving what is essentially a five-dimensional problem (three spatial dimensions and two angular dimensions). Improvements in computational resources have resulted in higher-resolution simulations for both climate and weather modeling, and this increase in resolution makes the assumption of purely vertical radiative transfer more and more difficult to justify. However, these improvements to computing power make it possible to attempt to solve a three-dimensional radiative transfer equation as part of a larger atmospheric model and in doing so take into account the three-dimensional effects that are commonly neglected. Doing this efficiently requires the use of adaptive mesh refinement, which while commonly used in the CFD community, is not a technique that is widely used in simulations of the atmosphere. In this work, Jexpresso.jl, an open source flexible general conservation law solver, is extended to solve the full equations of Radiative transfer. Adaptive mesh refinement in both space and angle is then used
to make it possible to speed up solving the three-dimensional radiative transfer equation and couple the solution to an atmospheric model which uses adaptive mesh refinement around clouds. The objective is to allow for improvements in the prediction of atmospheric flow behaviors around clouds and the resulting radiative heat fluxes, which are highly sensitive to cloud cover.
How to cite: Tissaoui, Y., Stechmann, S., Wang, H., and Marras, S.: Toward three-dimensional radiative transfer computations via adaptive mesh refinement in both space and angle, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-13120, https://doi.org/10.5194/egusphere-egu26-13120, 2026.
Various grids are being used to discretize the sphere for general circulation models of ocean and atmosphere, coupled to sea ice and land models. Depending on the numerics (or the network architecture for machine learning-based models) some grids are better suited and more widely used but no single grid meets all requirements. The HEALPix grid was originally invented for applications in cosmology and designed for spectral transform efficiency, hierarchical nested ordering and equal-area grid cells. Here, we present several variants: The OctaHEALPix grid, and the octaminimal Gaussian and Clenshaw-Curtis grids. The OctaHEALPix grid inherits all the properties of the original HEALPix grid but based on a single square matrix instead of 12 used for the base pixels (faces) of the HEALPix grid. This eliminates singularities for the 8 corner cells with 7 instead of 8 neighbouring cells. The nested hierarchy of the OctaHEALPix grid is therefore a quadtree from the coarsest to the finest zoom level. Data on the OctaHEALPix grid can be arranged as a single square matrix, representing an equal-area map projection of the sphere, allowing for interpolation-free visualisation and data storage. However, the HEALPix grids do not provide an exact quadrature in the Legendre transform, such that transform errors are higher than with Gaussian, or equi-angle latitude-based grids using the Clenshaw-Curtis quadrature. The inexact transform with the HEALPix grids does not pose any problems in simulations where other sources of error dominate. We present the grid-flexible spectral transform implemented in the atmospheric circulation model SpeedyWeather.jl that simultaneously supports all ring-based, equi-longitude grids both on CPU (including multithreading) and GPU. The OctaHEALPix grid is also favourable for machine learning-based models like diffusion models based on the UNet architecture which only require custom boundary conditions. The other HEALPix variant we present is the octaminimal Gaussian grid, which imposes Gaussian latitudes on the OctaHEALPix grid which reduces transform errors while preserving more of the HEALPix’s equal-area and hierarchy properties. Similarly, the octaminimal Clenshaw-Curtis grid uses regular latitudes with the Clenshaw-Curtis quadrature. Simulations based on these grids are presented for coupled climate simulations and implications for hybrid numerical and machine learning-based models are discussed.
How to cite: Klöwer, M., Gelbrecht, M., Leland, J., Groenke, B., and Hotta, D.: Variants of HEALPix grids for global climate modelling, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-15163, https://doi.org/10.5194/egusphere-egu26-15163, 2026.
We introduce a new suite of idealized test cases for the dynamical cores of atmospheric General Circulation Models. The tests were developed for the 2025 Dynamical Core Model Intercomparison Project (DCMIP-2025) which took place at the National Center for Atmospheric Research (NCAR) in June 2025 alongside a summer school (https://dcmip.org). The test suite was designed to probe (a) the impact of topography on atmospheric flows at various scales (see also the EGU 2026 companion paper led by Timothy Andrews); (b) a convection (squall line) test case; and (c) idealized experiments for atmospheric machine learning emulators.
This presentation focuses on the squall line test case, which challenges nonhydrostatic dynamical cores to accurately simulate moist convection and squall line dynamics at kilometer-scale resolutions. The scenario incorporates simplified moisture feedbacks using a warm-rain Kessler parameterization and is implemented across a range of horizontal (0.25–4.0 km) and vertical (250–500 m) grid spacings on a reduced-radius sphere. We assess three leading nonhydrostatic dynamical cores: the nonhydrostatic version of the Department of Energy/NCAR ‘Spectral Element’ model (called HOMME or SE), the ‘Model for Prediction Across Scales’ (MPAS), and NOAA’s Finite-Volume cubed-sphere dynamical core FV3. These are all available within NCAR’s Community Atmosphere Model (CAM) which is the atmospheric component of the Community Earth System Model (CESM).
Our analysis examines the numerical convergence characteristics, model-to-model variability, and the influence of core-specific dissipation mechanisms on the simulated convective storms. The results demonstrate the test case’s potential as a rigorous benchmark for future core development and model intercomparison efforts. This work contributes to advancing robust evaluation frameworks for atmospheric models in the kilometer-scale, nonhydrostatic regime.
How to cite: Jablonowski, C., Androski, N., Andrews, T., and Hughes, O.: DCMIP-2025: Introducing an Idealized Squall Line Test Case for Nonhydrostatic Dynamical Cores, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-8455, https://doi.org/10.5194/egusphere-egu26-8455, 2026.
Spectral methods have a well-established history in numerical weather prediction. However, stable formulations for high-order discontinuous Galerkin (DG) methods in both horizontal and vertical directions have only recently been developed using split-form DG (cf. Waruszewski et al. and Souza et al.).
This presentation provides an overview of the implementation of split-form DG methods for the compressible Euler equations, utilizing various formulations on both triangular and quadrangular spherical grids. Key focus areas include split-form versions on triangular grids and the linear algebra associated with HEVI-like (Horizontally Explicit, Vertically Implicit) temporal integration schemes. To manage vertically propagating sound waves implicitly, temporal integration is performed using W-Rosenbrock methods with an approximate Jacobian in the radial direction.
The framework is implemented in the Julia package CGDycore.jl, leveraging KernelAbstractions.jl and MPI for parallel execution across diverse architectures. We present numerical comparisons using the Held-Suarez test case, alongside comprehensive weak and strong scaling results for different computing environments.
How to cite: Knoth, O. and Artiano, M.: A discontinuous Galerkin weather dycore for triangular and quadrangular grids, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-14161, https://doi.org/10.5194/egusphere-egu26-14161, 2026.
Atmospheric dynamical cores solve the equations of motion that describe the fluid flow at the resolved scales. Commonly used dynamical cores, including the Department of Energy's (DOE) ‘Higher Order Methods Modeling Environment’ (HOMME) dynamical core, also known as a ‘Spectral Element’ (SE) dynamical core, and the National Center for Atmospheric Research's (NCAR’s) ‘Model for Prediction Across Scales’ (MPAS), contained two approximations of the fluid dynamics equations. They are called the Shallow-atmosphere and Traditional (SA+T) approximations. These approximations discard the so-called Nontraditional Coriolis terms. Omitting these terms significantly biases the dynamical response to diabatic heating, inducing up to 10% relative error in the wind field as linear model results suggested. The biases induced by these approximations are expected to be most severe in tropical regions, and to become more severe at finer grid spacings. We present preliminary evidence that removing these approximations may help mitigate the double Intertropical Convergence Zone (ITCZ) bias that is present in many climate models.
We present dynamical core intercomparisons between the shallow-atmosphere version of the HOMME and MPAS dynamical cores and their deep-atmosphere variants which remove the SA+T approximations. We present an overview of the modifications made to the HOMME and MPAS dynamical cores. Aquaplanet simulations made using NCAR’s Community Earth System Model (CESM) show that for sea surface temperature distributions that produce a double ITCZ, removing the SA+T approximations induces equatorward shifts in precipitation. We will examine sensitivities to physics-dynamics coupling strategies and parametric sensitivity to the deep convection scheme. We also examine the performance of MPAS, which collocates physics columns with dynamics columns, with HOMME, which simulates physics on a sparser finite-volume grid. Our simulations are done at a nominal 1º grid spacing, providing conclusive evidence that the SA+T approximations induce climatologically significant biases at the coarser grid spacings used in workhorse climate model configurations.
How to cite: Hughes, O., Ong, H., Herrington, A., Lauritzen, P., Guba, O., Taylor, M., and Jablonowski, C.: Examining the Contribution of Deep-Atmosphere Dynamical Cores to Mitigating Double ITCZ Bias in Aquaplanets, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-15662, https://doi.org/10.5194/egusphere-egu26-15662, 2026.
Normal modes of hydrostatic primitive equations on the sphere decompose the linearized flow into the slow propagating Rossby, fast propagating inertio-gravity (IG), and equatorial Kelvin and MRG waves, characterized by intermediate speeds. However, this simple characterization is no longer valid in the nonlinear system, where wave-wave, wave-mean flow interactions, and adiabatic forcing can significantly alter the wave speeds.
A number of methods aimed at determining the slowly evolving component of the nonlinear flow (so called slow manifold) were proposed under the umbrella terms ”nonlinear normal mode decomposition” (NNMD) or "nonlinear normal mode initialization" (NNMI). In the classical NNMD, Rossby waves are considered slow a priory, while the slow IG part consists of the the unbalanced component of the flow slaved to the Rossby modes. More precisely, in classical NNMD, slow IG waves are those that would be stationary in a nonlinear flow consisting of the slow modes only. While this approach is very successful in suppressing high-speed gravity waves, it also suppresses slowly propagating linearly unbalanced flow and large scale tropical circulation.
We propose a new method of flow decomposition into the slow and fast components based on computing instantaneous phase speeds of normal modes in the nonlinear system. The method does not make assumptions on the composition of the slow manifold. Instead, the decomposition is obtained from a constraint optimization problem that minimizes the norm of the fast component, while requiring that the slow manifold does not contain modes propagating faster than the selected cutoff speed. We apply the method to reanalysis data, demonstrate its efficiency, and analyze the composition of the slow manifold as function of the cutoff speed. In particular, when the cutoff is chosen to be approximately equal to the fastest linear Rossby wave speed, most of the tropical circulation and significant part of unbalanced modes are retained in the slow manifold.
How to cite: Vasylkevych, S., Chen, J., Holube, K., Zagar, N., and Lunkeit, F.: New nonlinear normal modes flow decomposition based on instantaneous phase speeds in the normal modes framework, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-15174, https://doi.org/10.5194/egusphere-egu26-15174, 2026.
Weather and climate models require atmospheric transport to be numerically stable, accurate, efficient, and mass-conserving. Adaptively Implicit-Explicit (AdImEx) time stepping can provide mass conservation and stability while taking long time steps for efficiency. However, challenges persist: (1) strongly reduced accuracy for large Courant numbers, and (2) a lack of preservation of uniform fields with substage time step sizes varying in space. We introduce AdHImEx, a novel AdImEx approach that resolves these issues, delivering an unconditionally stable, conservative, efficient advection scheme that preserves unity and increases temporal accuracy for large Courant numbers from first- to second-order. Remarkably, it requires only a single implicit solve per time step, performed with an iterative matrix solver, while retaining the full efficiency and third-order accuracy of the explicit time stepping in regions with small Courant numbers. Two-dimensional results will be shown.
How to cite: te Winkel, A., Weller, H., Kühnlein, C., and Kent, J.: AdHImEx: Adaptively High-order Implicit-Explicit Time Stepping for Atmospheric Transport with Long Time Steps, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-3905, https://doi.org/10.5194/egusphere-egu26-3905, 2026.
Accurate and efficient time integration is a critical component of Numerical Weather Prediction (NWP) and climate modelling. Current approaches in the Met Office’s next‑generation dynamical core, GungHo, rely on a serial-in-time, low-order, iterative semi‑implicit timestepping scheme coupled with Flux‑Form Semi‑Lagrangian (FFSL) transport scheme. As model resolutions and computational scales increase, exploiting parallelism in the time dimension is becoming increasingly important.
This work evaluates the potential of time‑parallel Deferred Correction (DC) methods as a viable alternative to these traditional schemes. I examine two strategies for introducing temporal parallelism:
- Revisionist Integral Deferred Correction (RIDC), which parallelises across correction sweeps, and
- Spectral Deferred Correction (SDC) with time‑parallel (diagonal) preconditioners, which parallelises across collocation nodes.
Results are presented for a suite of standard dynamical core test cases for the compressible Euler equations. These experiments demonstrate how time‑parallel DC algorithms can improve temporal accuracy and offer meaningful opportunities for reducing wall‑clock time.
Together, these developments highlight the promise of DC-based integrators as a pathway toward more scalable and efficient time-stepping in next‑generation atmospheric models.
How to cite: Brown, A.: A comparison of time-parallel deferred correction methods for atmospheric modelling, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-19200, https://doi.org/10.5194/egusphere-egu26-19200, 2026.
In the last decade, there has been growing interest in developing new dynamical cores for climate and weather simulations based on the Discontinuous Galerkin (DG) approach. To achieve accuracy, efficiency, and stability, various numerical formulations of the Euler equations are currently being explored.
In this talk, we present novel structure-preserving methods for different formulations of the Euler equations. These include models based on different thermodynamic variables, such as potential temperature, internal energy, or total energy, as well as the Exner pressure formulation and the vector-invariant form. These methods are developed within the flux-differencing DG framework, with a specific focus on the efficient implementation of split forms for both conservative and non-conservative terms.
The implementation is carried out entirely in Julia and is integrated within the Trixi.jl and CGDycore.jl ecosystems. We will discuss the different architectural approaches used in these packages and showcase numerical results, specifically focusing on the baroclinic instability test case to compare the different formulations.
How to cite: Artiano, M., Knoth, O., Spichtinger, P., and Ranocha, H.: Structure-Preserving Methods for the Euler Equations, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-18731, https://doi.org/10.5194/egusphere-egu26-18731, 2026.
The atmosphere, ocean, and sea ice components in Earth system models are coupled via boundary conditions at the sea surface. Standard coupling algorithms correspond to the first step of an iteration, so-called Schwarz waveform relaxation. Not iterating is computationally cheap but introduces a numerical coupling error, which we aim to quantify for the case of a coupled single column model: the EC-Earth AOSCM, which uses the same coupling setup and model physics as its host model, EC-Earth. To this end, we iterate until a reference solution is obtained and compare this with standard, non-iterative algorithms. Understanding the convergence behavior of the iteration, as well as the size of the coupling error, can inform model and algorithm development. Past studies demonstrated that the SWR solution eliminates phase errors, reduces ensemble spread, and can indicate whether current coupling setups are mathematically consistent. Our implementation is based on the OASIS3-MCT coupler and allows to estimate the coupling error of multi-day simulations.
In the absence of sea ice, SWR convergence is robust. Coupling errors for atmospheric variables can be substantial. When sea ice is present, results strongly depend on the model version. In the latest model version, coupling errors in sea ice surface and atmospheric boundary layer temperature are often large. Generally, we find that abrupt transitions between distinct physical regimes in certain parameterizations can lead to substantial coupling errors and even non-convergence of the iteration. We attribute discontinuities in the computation of atmospheric vertical turbulence and sea ice albedo as sources for these problems. We conclude the talk with some new theoretical results on analytically describing atmosphere-ocean-sea ice coupling.
How to cite: Schüller, V., Lemarié, F., Birken, P., and Blayo, E.: Quantifying Coupling Errors in Atmosphere-Ice-Ocean Coupling, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-287, https://doi.org/10.5194/egusphere-egu26-287, 2026.
In some geophysical flow phenomena, such as landslide-generated tsunamis and slow earthquake-generated waves, non-hydrostatic pressure has been shown to be crucial, resulting in an effect known as the dispersive effect. One practical approach to include such an effect is by extending the Shallow Water Equations (SWE), which can be achieved by splitting the pressure terms into hydrostatic and non-hydrostatic pressure while deriving a depth-averaged form. In the end, this model requires us to solve an elliptic system of equations to make a correction to the hydrostatic SWE approximation. However, solving the elliptic problem burdens computing time significantly. Therefore, we introduce a locally adaptive non-hydrostatic model that allows us to solve the extension locally. To achieve reliable results, the corrections need to be adapted in the area where the non-hydrostatic pressure might hold a significant role. We define these areas with a simple criterion based on the hydrostatic solution. To validate our model, we apply our model to a test case that involves moving bottom-generated waves. Our result shows that we can achieve a good agreement between the local and global models, where the former approach clearly reduces the computational time.
How to cite: Firdaus, K. and Behrens, J.: Two-Dimensional Locally Adaptive Non-Hydrostatic Model for Moving Bottom-Generated Waves, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-1213, https://doi.org/10.5194/egusphere-egu26-1213, 2026.
The Dynamical Core Intercomparison Project (DCMIP) is an endeavour to evaluate and compare dynamical cores through the application of idealised test cases. Previous editions of DCMIP, in 2008, 2012, and 2016, introduced many widely used cases, including baroclinic wave, tracer transport, supercell, and tropical cyclone tests. We present a new test case developed for the most recent edition of DCMIP, held in June 2025 (DCMIP2025). This case uses an easily implementable initial condition, with the addition of mountain orography, to generate nonlinear mesoscale flow phenomena. Two different orographies are considered, which lead to dynamics of a gap flow and vortex shedding. The gap flow case exhibits accelerated flow and the creation of a pair of symmetric lee vortices, whilst the vortex shedding case instigates a flow reminiscent of a von Kármán vortex street. These tests are compared in four dynamical cores: three from the United States, along with GungHo, the compatible finite element dynamical core from the UK Met Office’s new LFRic model. Contrasting simulations with each dynamical core highlights differences in the numerical designs, including the grid choice and sources of numerical diffusion.
How to cite: Andrews, T., Jablonowski, C., Hughes, O., and Bendall, T.: New mountain-generated mesoscale flow test cases from DCMIP2025, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-2009, https://doi.org/10.5194/egusphere-egu26-2009, 2026.
Rayleigh-Bénard convection (RBC) — the buoyancy-driven instability that arises when a fluid is heated from below and cooled from above, leading to the formation of convective plumes — serves as a fundamental and challenging benchmark problem in geophysical fluid dynamics. Simulating RBC at high Rayleigh numbers demands extremely fine spatial resolution, which in turn requires high-performance computing (HPC) resources to achieve results within feasible runtimes. However, strong parallel scaling based on spatial domain decomposition eventually saturates due to communication overheads, and simulations involving very large numbers of time steps can still be prohibitively time-consuming, with limited scope for further speedup through additional spatial parallelism.
To overcome these limitations, time-parallel algorithms — which introduce concurrency along the temporal dimension — offer a promising approach to extend strong scaling beyond the saturation of space-only parallelization. Despite their potential, constraints imposed by causality often make these methods challenging to design, and some classes of parallel-in-time algorithms can suffer from poor parallel efficiency. In contrast, parallel-across-the-method techniques, while providing only small-scale parallelism, tend to be easier to implement and can achieve competitive efficiency. Previous efforts to develop parallel Runge-Kutta methods were only moderately successful, primarily because the associated stability restrictions were more stringent than those of their serial counterparts.
Parallel Spectral Deferred Corrections (pSDC), however, enable parallel computation of stages without significantly reducing the maximum stable time step. In this talk, we introduce pSDC and present a bespoke solver that combines pSDC with parallel Fast Fourier Transforms (FFTs) on GPUs to facilitate efficient, large-scale simulations of RBC. We show performance results obtained on the JUWELS Booster and JUPITER HPC systems at the Jülich Supercomputing Centre, showcasing how pSDC can achieve runtime reductions beyond the limits of spatial parallelization alone.
How to cite: Ruprecht, D., Saupe, T., Lunet, T., Götschel, S., and Speck, R.: Improved parallel scaling of 3D Rayleigh-Benard convection simulations by parallelization in time, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-3094, https://doi.org/10.5194/egusphere-egu26-3094, 2026.
We propose a new formulation of a three-dimensional spectral model for the primitive equations, in which the spectral method is applied not only in the horizontal but also in the vertical direction. In this model, we utilize scaled Laguerre functions as the vertical basis and spherical harmonics as the horizontal basis. This formulation introduces a tunable scaling parameter, enabling the model top to be placed at high altitudes while maintaining adequate vertical resolution in the upper atmosphere. This formulation is implemented as a numerical model, and its performance is validated through a series of benchmark experiments. The results demonstrate that the numerical error of the present model decreases much faster than that of a corresponding vertical finite-difference model as the vertical degrees of freedom are increased. Furthermore, using linearized two-dimensional primitive equations, the properties of gravity waves under the proposed discretization are examined and compared with those of conventional vertical finite-difference discretizations. The results confirm that the proposed method can accurately represent upward-propagating waves.
How to cite: Fujita, S. and ishioka, K.: Development of a Three-dimensional Spectral Atmospheric Dynamical Core using Scaled Laguerre Functions, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-3274, https://doi.org/10.5194/egusphere-egu26-3274, 2026.
The choice of grid staggering has been a fundamental design decision in atmospheric dynamical core development for decades. Conventional wisdom, largely derived from low-order Von Neumann analysis, holds that the unstaggered A-Grid exhibits poor dispersion properties near the grid scale, making it unsuitable for geophysical fluid dynamics. This perception has steered the community toward staggered grid formulations (B-, C-, or D-Grid) despite the algorithmic complexity they introduce.
We present a numerical approach to Von Neumann analysis that enables rigorous evaluation of high-order schemes with complex time-stepping methods, including implicit and semi-implicit formulations such as the forward-backward scheme. By numerically solving the Fourier-transformed equations across the two-dimensional space of Courant numbers and numerical phase, this method circumvents the algebraic intractability that has limited traditional analysis to simplified low-order cases.
Our analysis reveals a critical finding: the dispersion and dissipation differences attributable to grid staggering choices diminish substantially with high-order spatial discretization. When combined with the Low Mach number Approximate Riemann Solver (LMARS), which provides implicit scale-selective diffusion matched to the dispersion characteristics, the A-Grid formulation effectively controls numerical noise while maintaining accuracy for well-resolved modes. Idealized tests with discontinuous wave packets validate these theoretical predictions and demonstrate that high-order LMARS produces significantly less numerical noise than inviscid schemes on both staggered and unstaggered grids.
These findings carry significant implications for next-generation dynamical core design. The A-Grid formulation offers compelling advantages: algorithmic simplicity facilitating GPU implementation, straightforward conservation properties, collocated variables simplifying physics-dynamics coupling, and natural compatibility with data-driven approaches in hybrid modeling. Continued adherence to conventional wisdom rooted in low-order analysis risks misguiding the development of dynamical cores optimized for modern computing architectures and emerging AI-integrated Earth system models.
How to cite: Chen, X.: Rethinking A-Grid for next-generation dynamical cores: High-order numerical analysis challenges conventional wisdom, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-4623, https://doi.org/10.5194/egusphere-egu26-4623, 2026.
The multiscale finite element method (MsFEM) was originally proposed for stationary (elliptic) or quasi-stationary (parabolic) problems. It was extended to linear transport-dominated problems, utilizing a semi-Lagrangian subgrid reconstruction (SLMsR) approach. In this presentation we introduce the extension of the method to coupling non-linear systems of hyperbolic equations. To demonstrate the accuracy and applicability of SLMrS we use a shallow water model, where we discretize the momentum equation on a fine mesh and inform the coarse-mesh continuity equation of subgrid-scale features by using a multiscale basis. We show accuracy and convergence of the new method for smooth and non-smooth test cases.
How to cite: Behrens, J. and Gabsi, M.: Multiscale finite element method in a shallow water model, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-4907, https://doi.org/10.5194/egusphere-egu26-4907, 2026.
We investigate the coupling of nonlinear water-wave motion to heaving buoy wave-energy dynamics [1] in the presence of an inequality constraint. Building on augmented Lagrangian variational principles (VPs) developed by Burman [2] and others, we impose constraints of the form G(q)≥0, where q involves system variables, through a Lagrange multiplier λ. The strict Kuhn–Karush–Tucker (KKT) conditions {λ G=0,G(q)≥0, λ≤0} are replaced by those smooth approximations of the involved function F(c G(q)−λ)=max(c G(q)−λ,0), with smoothing parament c>0, allowing explicit computation of the multiplier λ as (part of) a force. Our approach combines: (a) an Average Vector Field (AVF) energy-conserving time-stepping method, extended to water-wave systems with an auxiliary field, enforcing energy conservation in the discrete system; (b) a (novel) smooth relation λ(G) that regularises the KKT conditions by approximating the solution G=0 with λ≤0 and G>0 with λ=0 in the (λ,G)-plane, but leading to an implicit definition of the function F(c G(q)−λ). This framework has been implemented and tested in the finite-element environment Firedrake, leading to improved and surviving benchmarks for the problems: (i) a point particle under gravity bouncing off a rigid table, (ii) a particle moving in a rectangular (“billiard”) domain, and (iii) forced (Variational “Boussinesq” Model-type) nonlinear water waves in a horizontal channel causing buoy motion in a wave-enhancing contraction. The latter, finite-element, model supports design of a prototype wave-energy device for enhanced energy capture. More generally, this work aims to develop analytical and computational tools for finite-element coupling of nonlinear wave dynamics in fluid-structure interactions, here exemplified by the vertical (heave) motion of the buoy.
[1] O. Bokhove, A. Kalogirou and W. Zweers (2019) From Bore–Soliton–Splash to a new wave-to-wire wave-energy model. Water Waves 1, 217-258.
[2] E. Burman, P. Hansbo and M.G. Larson (2023) The augmented Lagrangian method as a framework for stabilised methods in computational mechanics. Archives of Computational Methods in Eng. 30, 2579–2604.
How to cite: Bokhove, O.: Fluid-structure interactions with numerics for a smoothed augmented Lagrangian variational principle applied to a wave-energy device, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-10650, https://doi.org/10.5194/egusphere-egu26-10650, 2026.
Two-fluid magnetohydrodynamics (MHD) extends single-fluid MHD by retaining separate ion and electron thermodynamics together with an extended Ohm’s law. This enables the capture of Hall physics, dispersive waves, and fast magnetic reconnection – phenomena that are inaccessible to single-fluid MHD and are central to applications such as magnetic confinement fusion, the heliosphere, and the Earth’s magnetosphere.
In this presentation, we discuss a novel spatial discretization for the two-fluid MHD equations based on compatible finite element methods. The approach preserves important structural properties including the divergence-free constraint on the magnetic field, energy conservation, and a consistent treatment of both fluid and magnetic helicity. In particular, the ion velocity and magnetic field spaces are chosen to admit a natural discrete definition of the diagnostically determined electron velocity.
The resulting scheme is designed for low-dissipation regimes in which fluid transport dominates, a setting relevant to many scenarios of interest in the aforementioned applications. To ensure a stable field evolution, we incorporate transport stabilization for all fields while preserving the underlying structural properties of the two-fluid MHD system. The stabilization is based on interior penalty methods, extending our previous work on magnetic field transport in single-fluid resistive MHD (Wimmer, Tang, 2024). Numerical experiments demonstrate the structure preserving and stabilization properties of the method through test cases focusing on fluid helicity, magnetic helicity, and the relative decay rates of helicity and energy in the presence of dissipation.
References
Golo A. Wimmer and Xian-Zhu Tang (2024), Structure preserving transport stabilized compatible finite element methods for magnetohydrodynamics, Journal of Computational Physics, Volume 501, 112777.
How to cite: Bauer, W., Wimmer, G. A., and Tang, X.-Z.: Structure preserving, transport stabilized finite element methods for two-fluid magnetohydrodynamics, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-13095, https://doi.org/10.5194/egusphere-egu26-13095, 2026.
Data-driven algorithms provide an exciting opportunity to represent unresolved, and therefore parameterised, processes in a way that matches available data without the need to `hand-tune' the parameterisation. However, on their own they suffer from issues with accurate prediction of extreme events and long-term trends, at least in part due to their lack of any representation of physical constraints. Hybrid weather and climate prediction models address this issue by combining data-driven algorithms with more traditional algorithms based on solving the partial differential equations (PDEs) that govern atmospheric flow. However, training the data-driven component separately from the PDE solver can introduce issues with stability and still does not ensure that the combined model will not drift over long times. Online training addresses this issue by more tightly coupling the two model components during training. This requires that the PDE model is differentiable so that during training backpropagation can be performed on multiple timesteps of both the PDE solver and the data-driven component. In this work we introduce a compatible finite element dynamical core that is automatically differentiable since it is built using the Firedrake finite element library. Here we will show how this can be coupled to PyTorch to perform end-to-end training on idealised test cases.
How to cite: Shipton, J. and Hartney, N.: Investigating the impact of model formulation on the accuracy andefficiency of a hybrid dynamical core with online training., EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-20183, https://doi.org/10.5194/egusphere-egu26-20183, 2026.
We present a novel staggered discontinuous Galerkin method on general polytopal meshes for solving the shallow water equation on the sphere. Every polygon is decomposed into so called kite quads by joining the vertices with the midpoints of the adjacent edges and the midpoint of the polygon. The shallow water equation is solved in vector-invariant form whereby vorticity is determined diagnostically. Height and momentum are discretized on the primal respectively dual cells whereby the composed finite elements are continuous over internal edges and discontinuous over boundary edges. This results in block-diagonal mass matrices. Mass lumping can reduce the fill in of these matrices further.
The method is implemented in Julia in the package CGDycore.jl which unifies different numerical dycores under one umbrella. Numerical results are presented for standard test cases on different spherical grids.
How to cite: Knoth, O.: Staggered finite element methods on polytopal meshes, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-10749, https://doi.org/10.5194/egusphere-egu26-10749, 2026.
