Validation of an OpenFOAM^®-based solver for the Euler equations with benchmarks for mesoscale atmospheric modeling

Forecasting the rapid changes in the Earth’s climate is one of the biggest challenges of our times. Since the complexity of the problem requires inputs and close collaboration from scientists across various disciplines, ranging from physics of the atmosphere to computer science, open-source software has emerged as the evident choice. In fact, the availability of open-source software for use, modification, and distribution makes it ideal for collaborative development. Some examples of open-source software for climate simulations are the Community Earth System Model,¹ the Energy Exascale Earth System Model,² the MIT General Circulation Model,³ and the Climate Machine.⁴ All of these examples have large communities of active users and developers that have contributed over several years. Since climate simulations are an extension of weather forecasting, we note that open-source software has been a tool of choice for weather prediction as well. Examples include WRF⁵ and Atlas.⁶ 

Fast and accurate weather/climate forecasts need state-of-the-art numerical and computational methodologies. Although the above-mentioned software packages are great tools for realistic simulations, testing and assessing new numerical approaches within them are non-trivial. In fact, the very complexity that allows for realistic simulations requires a considerable amount of time to familiarize with, making such software impractical as a testbed. This paper lays the foundation for a new open-source package, called Geophysical and Environmental Applications (GEA),⁷ specifically created for the assessment of new computational approaches for the simulation of mesoscale atmospheric flows and ocean flows.^8,9 This idea is that if a new approach fails to meet desired accuracy or efficiency standards in our simplified software, it would not be considered for implementation in more advanced software. To maximize the reach and impact, we choose to build our software package on OpenFOAM^®,¹⁰ an open-source, freely available C++ finite volume library that has become widely used in Computational Fluid Dynamics (CFD).

In recent years, OpenFOAM has been adopted to study a large variety of computational approaches for different CFD applications, including plasma cutting,¹¹ fire plumes,¹² and metal forming processes.¹³ Thanks to parallel computing support, multiphase modeling capabilities, a wide range of existing turbulence models, and ease of implementation for new ones, OpenFOAM represents a great tool for the development and assessment of computational strategies for atmospheric modeling. In this context, researchers have chosen OpenFOAM to perform atmospheric boundary layer simulations (see, e.g., Refs. 14–16), multi-fluid atmospheric convection (see, e.g., Ref. 17), and mesh r-adaptivity for weather prediction over steep terrain (see, e.g., Ref. 18). The results in Refs. 17 and 18 have been obtained with AtmosFOAM,¹⁹ an OpenFOAM-based model of the global atmosphere using arbitrary-shaped cells. Our solver starts with a focus on local atmospheric problem. However, we do not exclude interfacing it with AtmosFOAM in the future.

As the core of our open-source package, we present an OpenFOAM-based solver for the Euler equations for stratified fluid flow, which are of importance in non-hydrostatic mesoscale atmospheric modeling, and assess it through two well-known test cases. We consider the Euler equations written in conservative form using density, momentum, and total energy as variables. It is shown in Ref. 20 that this set of equations yields less dissipative results than other forms of the Euler equations for relevant test problems. Hence, it was recommended to pursue this formulation for the development of mesoscale models, although other formulations (i.e., non-conservative forms and the conservative form using density, momentum, and potential temperature as variables) have received more attention for atmospheric studies. For the stabilization of the Euler equations and to capture sub-grid processes, we will assess and compare two Large Eddy Simulation (LES) models built in OpenFOAM: the classical Smagorinsky model²¹ and the one equation eddy-viscosity model.²² We will show that both LES models yield numerical results that agree very well with data published in the literature.

Although, currently, both pressure-based and density-based approaches for the Euler equations apply to a wide range of flows, from incompressible to highly compressible, pressure-based solvers have been conventionally used for mildly compressible flows (as is the case for atmospheric flows). Despite this, the majority of the software for weather prediction adopts density-based approaches, which were originally designed for high-speed compressible flows. In this paper, we opt for a pressure-based solver and show that it yields accurate results when compared to data in the literature. In the spirit of OpenFOAM, our solver uses a splitting scheme that decouples the computation of each variable in order to achieve high computational efficiency.

The rest of this paper is organized as follows: In Sec. II, we briefly describe the formulation of the Euler equations and the LES models under consideration. Section III presents our pressure-based approach and provides the details of space and time discretization. Numerical results for the two benchmark tests are discussed in Sec. IV. Conclusions are drawn in Sec. V.

We consider the dynamics of dry atmosphere in a spatial domain of interest Ω by neglecting the effects of moisture, solar radiation, and heat flux from the ground. We assume that dry air behaves like an ideal gas.

The numerical discretization of the model presented in Subsection II A would lead to an accurate description of atmospheric flow if one could afford a discretization mesh able to capture all the scales of turbulent structures. Since the Kolmogorov scale of a typical atmospheric problem is about 10⁻⁴ m, a Direct Numerical Simulation is currently beyond reach. One way to keep the computational cost affordable is to adopt a coarser mesh and model the effects of the small scales that are not directly solved. This can be carried out via Large Eddy Simulation.

This section presents space and time discretization for the generic LES models (1), (4), (5), (7), (11), and (12).

A monolithic approach for coupled problems (15)–(21) would lead to high computational costs. Thus, to save computational time, we adopt a splitting approach. While the splitting applies to the partial differential equations, non-differential equations in problems (15)–(21) are used as predictors or to compute end-of-step variables. The splitting scheme is as follows: given ρ⁰, u⁰, h⁰, p⁰, and T⁰, set K⁰ = |u⁰|²/2 and for n ≥ 0 perform

• Step 1: First, find intermediate density $ρn+13$⁠, intermediate velocity $un+13$⁠, and associated kinetic energy density $Kn+13$ such that

  $ρn+13Δt=bρn+1−∇⋅(ρnun),$

  (22)

  $ρn+13un+13Δt+∇⋅ρnun⊗un+13+∇pn−∇⋅2μanϵun+13+∇23μan∇⋅un=bun+1,$

  (23)

  $Kn+13=|un+13|22.$

  (24)

  Note that Eq. (24) is a predictor step, and the associated corrector step is Eq. (32).
• Step 2: Find specific enthalpy hⁿ⁺¹, temperature Tⁿ⁺¹, and second intermediate density $ρn+23$ such that

  $ρn+13hn+1Δt+∇⋅(ρnunhn+1)−∇⋅μanPr∇hn+1=b̃en−ρn+13Kn+13Δt−∇⋅ρnunKn+13+pnΔt−ρn+13gun+13⋅k̂,$

  (25)

  $hn+1−cpTn+1=hn−cpTn,$

  (26)

  $ρn+23RTn+1=pn,$

  (27)

  where $b̃en=(ρnhn+ρnKn−1−pn−1)/Δt$⁠. Equation (26) will be used to compute the end-of-step temperature (see Sec. III A), whereas Eq. (27) is a predictor step for density. The associated correction will come from Eq. (30) as explained in Sec. III A.
• Step 3: Find end-of-step velocity uⁿ⁺¹ and associated kinetic energy density Kⁿ⁺¹, pressure pⁿ⁺¹ and pressure fluctuation p′^,n+1, and end-of-step density ρⁿ⁺¹ such that

  $ρn+13un+1Δt+∇⋅(ρnun⊗un+1)+∇pn+1−∇⋅(2μanϵ(un+1))+∇23μan∇⋅un=bun+1,$

  (28)

  $pn+1−p′,n+1=ρn+23gz,$

  (29)

  $pn+1−ρn+1RTn+1=0,$

  (30)

  $ρn+23Δt+∇⋅ρn+23un+1=bρn+1,$

  (31)

  $Kn+1=|un+1|22.$

  (32)

  Note that in Eq. (28), one diffusive term is kept implicit, whereas the other is treated explicitly. This is customary in OpenFOAM.

Section III A will describe the space discretization of the equations at each step and explain how to further reduce the computational cost associated with Step 3.

For the space discretization of Eqs. (22)–(32), we adopt a finite volume method. This requires partitioning the computational domain Ω into cells or control volumes Ω_i, with i = 1, …, N_c, where N_c is the total number of cells in the mesh. Let A_j be the surface vector of each face of the control volume, with j = 1, …, M.

We test our OpenFOAM-based solver for the Euler equation with two standard benchmarks that have been widely used to assess atmospheric dynamical cores, i.e., the rising thermal bubble and the density current. Both test cases involve a perturbation of a neutrally stratified atmosphere with uniform background potential temperature over a flat terrain. Therefore, before reporting the results for the two benchmarks, in Sec. IV A, we show that an unperturbed stratified atmosphere with uniform background potential temperature over a flat terrain remains unchanged up to a certain tolerance. In Sec. IV B, we present our results for the rising thermal bubble benchmark. There exist several variations of this benchmark, featuring different geometries and/or initial conditions. We use the settings from Ref. 30. See also Ref. 31 for a recent study using this variation. Our results for the classical density current test^32,33 are shown in Sec. IV C. Both benchmarks are two-dimensional. Although our solver is designed to deal with fully three-dimensional benchmarks, two-dimensional problems can be solved by generating a three-dimensional mesh with just one cell in the third dimension.

We would like to point out that neither the rising bubble nor the density current benchmark has an exact solution. Hence, one can only have a relative comparison with other numerical data available in the literature.

The goal of this first test case is to verify that an initial resting atmosphere over a flat terrain remains still within a reasonable accuracy for a long time interval.

The computational domain in the xz-plane is Ω = [0, 16 000] × [0, 800] m². In this domain, the hydrostatic atmosphere, initially at rest, is free to evolve until t = 25 days.³⁴ Impenetrable, free-slip boundary conditions are imposed at all the boundaries. We consider a uniform mesh with mesh sizes h = Δx = Δz = 250 m and we set the time step to Δt = 0.1 s.

Figure 1 shows the time evolution of the maximal vertical velocity w_max. We observe that w_max has an initial, rather fast, growth but subsequently stabilizes around 1 × 10⁻⁵ m/s. Thus, we conclude that the hydrostatic equilibrium is preserved with reasonably good accuracy, especially in view of the two benchmarks we are going to study next.

This benchmark is characterized by Re ≈ 1.8 × 10⁹ and M ≈ 0.04 when we set the characteristic quantities as follows: L_r = 2000 m (i.e., the initial bubble radius), ρ_r = 1.2 Kg/m³ and μ_r = 1.8 × 10⁻⁵ Pa s (i.e., the air density and viscosity at 300 K), T_r = 300 K, and U_r = 13.3 m/s (i.e., the maximum vertical component of the velocity).

For this test, we use five different meshes with uniform resolution h = Δx = Δz = 250, 125, 62.5, 31.25, 15.625 m. The time step is set to Δt = 0.1 s for all the simulations. For stabilization, we consider two strategies. First, following Ref. 30, we set μ_a = 15 and Pr = 1. Hereinafter, we refer to this model as AV15, where AV stands for artificial viscosity. We note that both of these are ad hoc values chosen by the authors of Ref. 30 to stabilize the numerical simulations. It is not unusual in benchmarks to set Pr = 1 although the air Prandtl number is about 0.71 at 20 °C (see, e.g., Ref. 36). Other authors have chosen other arbitrary values, like Pr = 0.1 in Ref. 37. For a qualitative analysis of the results for the rising thermal bubble as Pr varies we refer to the reader to Ref. 38. After AV15, we consider the Smagorinsky model as described in Sec. II B.

Figure 2 reports the perturbation of potential temperature θ′ at t = 1020 s computed by the AV15 model with all the meshes under consideration. By t = 1020, the air warmer than the ambient has risen due to buoyancy and deformed due to shearing motion. As a result, the bubble has evolved into a mushroom shape. From Fig. 2, we observe that the no substantial change in the computed θ′ when the mesh is refined past h = 125 m. We remark that to facilitate the comparison of the panels in Fig. 2, we have forced the color bar to range from 0 to 1. Qualitatively, these results are in very good agreement with those reported in the literature. For example, see Refs. 30, 35, and 38.

For further qualitative assessment, Fig. 3 displays velocity components u and w at t = 1020 s computed by the AV15 model with mesh h = 125 m. These contour plots are in very good agreement with those reported in the literature. See, e.g., Fig. 7 in Ref. 30.

For a more quantitative comparison, in Fig. 4, we show the time evolution of the maximum perturbation of potential temperature $θmax′$ and maximum vertical component of the velocity w_max computed by the AV15 model with all the meshes, together with the corresponding results from Ref. 30. We observe that the evolution of $θmax′$ computed with meshes h = 250 m and h = 125 m is affected by spurious oscillations. Such oscillation disappears at higher resolution and they are not present at all in the evolution of w_max. We see that $θmax′$ and w_max computed with meshes h = 31.25 m and h = 15.625 m are practically overlapped for the entire time interval, indicating that we are close to convergence. The “converged” w_max overlaps with the reference value until about t = 500 s and it remains close to it until about t = 800 s. The agreement of the “converged” $θmax′$ with the results from Ref. 30 is not as good. These trends are confirmed from Table I, which reports the extrema for the vertical velocity w and potential temperature perturbation θ′ at t = 1020 s obtained with the AV15 model, together with the values extracted from the figures in Ref. 30. The results from Ref. 30 are obtained with mesh resolution of 125 m and a density-based approach developed from a Godunov-type scheme that employs flux-based wave decompositions for the solution of Riemann problem. The authors of Ref. 30 start from the Euler equations written in density, velocity, and potential temperature. Other differences with our methodology include the orders of space and time discretization and the different treatment of the hydrostatic term. Given all these differences, we believe that our results cannot be considered off with respect to the reference. In addition, other references, such as Ref. 38 (see Fig. 1, top-left panel), report a $θmax′$ at t = 1020 s computed with resolution 125 m closer to 1 K (like in our case) than to 1.4 K (like in Ref. 30).

Finally, we run this test with the Smagorinsky model. We set the parameter C_s = 0.094, with C_k = 0.21 and C_ϵ = 1.048. Details on the choice of these parameters are reported in Sec. IV C. Figure 5 depicts the spatial distribution of θ′ at t = 1020 s computed with meshes 15.625 and 32.25 m. As expected, with these fine meshes, the Smagorinsky model can capture a larger amount of vortical structures than the AV15 model. The beautiful vortex strip created by the Rayleigh–Taylor instability at the edge of the bubble shown in Fig. 5 (left) matches well with that shown in Ref. 38 (see Fig. 1, bottom-right panel). Table II reports the extrema for the vertical velocity w and potential temperature perturbation θ′ at t = 1020 s obtained with the Smagorinsky model. Since we could not find any data obtained with an LES model for the exact setting of Ref. 30, Table II does contain reference data for comparison.

We conclude this section with a comment on the computational time required by the simulations to indicate the solver efficiency. We ran all the simulations using an 11th Gen Intel(R) Core(TM) i7-11700 @ 2.50 GHz processor. The full simulation with the coarsest mesh we considered (h = 250 m) takes about 14 s. It was only run with the AV15 model. With the finest mesh (h = 15.625 m), the full simulation with the AV15 model takes about 35 min, while it takes roughly 57 min with the Smagorinsky model.

This benchmark is characterized by Re ≈ 9 × 10⁹ and M ≈ 0.097 when we set the characteristic quantities as follows: L_r = 4000 m (i.e., the semi-major axis of the initial perturbation), ρ_r = 1.2 Kg/m³ and μ_r = 1.8 × 10⁻⁵ Pa s (i.e., the air density and viscosity at 300 K), T_r = 300 K, and U_r = 33.7 m/s (i.e., the maximum velocity magnitude at t = 900 s, which is the time the flow is fully developed).

We first consider uniform meshes with commonly used^{20,30,33,38,39} mesh sizes h = Δx = Δz = 400, 200, 100, 50, 25 m. The time step is set to Δt = 0.1 s. Since this benchmark features more complex vortical structures than the rising thermal bubble, we consider three stabilization strategies. Following Refs. 30 and 33, we first set μ_a = 75 in Eqs. (11) and (12) and Pr = 1 in Eq. (12). We will refer to this model as AV75 since 75 is an ad hoc artificial value. Then, we consider the Smagorinsky and kEqn models as described in Sec. II B.

Let us start with the AV75 model and a qualitative illustration of the flow evolution. Figure 7 shows θ′ computed with the AV75 model and mesh h = 25 m (i.e., the finest mesh among those considered) at t = 300, 600, 750, 900 s. We observe very good agreement with the results reported in Fig. 1 of Ref. 33, which were obtained with the same resolution. Indeed, we see that, as expected, the cold air descends due to negative buoyancy and strong downdrafts develop at the center of the cold bubble. When the cold air reaches the ground, it rolls up and forms a front. As this front propagates, shear is generated at its top boundary with a resulting Kelvin–Helmholtz-type instability that leads to a three-rotor structure at t = 900 for the resolution under consideration.

In Fig. 8, we report θ′ computed at t = 900 s by the AV75 model with all the meshes mentioned above. We observe the emergence of a clear three-rotor structure when the resolution is equal to or smaller than h = 100 m, with more definition when the mesh is finer. Also, the results in Fig. 8 are in very good agreement with those reported in the literature. See, e.g., Refs. 30, 33, ³⁸ and ³⁹.

Next, we focus on the Smagorinsky model. Given the poor quality of the solutions computed with the AV75 model and meshes coarser than h = 100 m, we will use only the meshes finer than h = 100 m. Moreover, following Ref. 38 for the LES models, we will also consider a very fine mesh with size h = 12.5 m. Figures 9, 10, and 11 display the time evolution of the potential temperature fluctuation computed with the Smagorinsky model and meshes h = 12.5, 25, 50 m, respectively. For all these simulations, the parameter C_s in Eq. (13) was set to 0.454, with C_k = 0.6 and C_ϵ = 1.048. We note that the ratio of the values of C_s used for this benchmark and for the rising thermal bubble is about 5, which is the ratio of the ad-hoc artificial viscosities. From Figs. 9–11, we clearly see that more vortical structures appear as the mesh size decreases. Indeed, we can see that with h = 50 m, there is a three-rotor structure at t = 900 s (see Fig. 11, bottom-right panel), and with h = 25 m, a quadri-rotor structure is present (see Fig. 10, bottom-right panel). The largest recirculation in the bottom-right panels of 11 and 10 is broken into two recirculations when mesh h = 12.5 m is used. See Fig. 9, bottom-right panel.

We now switch to the kEqn model. Figures 13, 14, and 15 display the time evolution of the potential temperature fluctuation computed with the kEqn model and meshes h = 12.5, 25, and 50 m, respectively. We observe a small difference between the solutions computed by the Smagorinsky model and the kEqn model with mesh h = 50 m for the entire duration of the time interval. Compare Fig. 11 with Fig. 15. The differences in the solutions given by the two LES models with mesh h = 25 m become visible around t = 750 s when the largest recirculation given by the kEqn model has a more flattened top boundary. In addition, the quadri-rotor structure at t = 900 s given by the Smagorinsky model is more defined. Compare the bottom panels in Fig. 10 with Fig. 14. This is expected since the kEqn model can be inaccurate for flows with strong streamline curvature and flows with large-scale turbulence structures as mentioned in Sec. II B. The differences in the solutions given the two LES models with mesh h = 12.5 m are pretty remarkable already at t = 600 s. Compare Fig. 9 with Fig. 13.

Figure 16 shows the time evolution of the space-averaged eddy viscosity (53) given by the kEqn model for meshes h = 12.5, 25, and 50 m. By comparing Fig. 16 with Fig. 12, we observe that the two LES models have a very similar evolution of μ_av for mesh h = 12.5 m. For the other two meshes, the Smagorinsky model ramps up the value of μ_av faster than the kEqn model. Finally, we note that μ_av introduced by both models keeps increasing over the time interval [0, 900] s since more and more vortical structures develop.

For a more quantitative comparison, we consider the front location, defined as the location on the ground where θ′ = −1 K, at t = 900 s. In Table III, we report the space interval that contains the front location for our implementations of the AV75, Smagorinsky, and kEqn models. We compare our values with the data in Ref. 33 (Table IV), which refer to the AV75 model and 14 different approaches for the numerical solution of this benchmark problem. We note that our results fall well within the values reported in Ref. 33. Furthermore, from Table III, we see that with our implementation of the Smagorinsky and kEqn models, the front becomes faster as the mesh is refined. This is more evident with the Smagorinsky model than with the kEqn model, while the AV75 model shows the opposite trend. The front locations obtained with the AV75 model and different meshes are within about 500 m from each other. The same is true for the kEqn model, while the Smagorinsky model gives a larger variation in front location (about 900 m) as the mesh size varies. Since, in Ref. 33, no results are reported for mesh h = 12.5 m, in Table III, we list the front location found with an LES model and resolution 12.5 m in Ref. 38. Such front location is within about 500 m of the value we obtain with the kEqn model and about 900 m of the value we obtain with the Smagorinsky model.

For a further quantitative comparison, we consider the potential temperature perturbation θ′ at t = 900 s along the horizontal direction at a height of z = 1200 m. Figure 17 displays a comparison between the results given the AV75 model on meshes h = 100, 50, 25 m and the results obtained with the same model and resolution 25 m in Ref. 20. Two numerical methods were used in Ref. 20 to approximate the solution: a spectral element method and a discontinuous Galerkin method, which provide curves that are superimposed in Fig. 17 and are denoted as “Reference.” We observe that our results are slightly out of phase for meshes h = 50, 100 m with respect to the data in Ref. 20. The potential temperature fluctuation obtained with mesh h = 25 m is in phase but has larger negative peaks, especially in correspondence of the first recirculation. However, given all the differences between our approach and the one in Ref. 20, the comparison is satisfactory.

As for the previous benchmark, we conclude by providing the computational times. The full simulation with the coarsest mesh (h = 400 m) takes about 14 s. It was run only with the AV75 model. The full simulation with the finest mesh (h = 25 m) takes about 53 min with the AV15 model, 59 min with the Smagorinsky model, and 1 h with the kEqn model.

The goal of this paper is to show that OpenFOAM provides a reliable and accurate framework for the simulation of non-hydrostatic atmospheric flows. To achieve this goal, we developed a pressure-based solver for the Euler equations written in conservative form using density, momentum, and total energy as variables and tested it against numerical data available in the literature for two well-known benchmarks, namely, the rising thermal bubble and the density current. For the rising bubble and the density current tests, we obtain good qualitative and quantitative comparisons when the classical Smagorinsky model and/or the one equation eddy-viscosity model are adopted for stabilization. We observed a small difference between the solutions computed by the Smagorinsky model and the one equation eddy-viscosity model for coarser meshes. The differences become visible for finer meshes, and the finer the mesh is, the earlier they appear. For both LES models, the front in the density current benchmark becomes faster as the mesh is refined. However, while the front locations obtained with the one equation eddy-viscosity model and different meshes are within about 500 m from each other, a larger variation in front location (about 900 m) is given by the Smagorinsky model as the mesh size varies.

We thank Dr. S. Marras for the fruitful conversations. We acknowledge the support provided by the European Research Council Executive Agency by the Consolidator Grant project AROMA-CFD “Advanced Reduced Order Methods with Applications in Computational Fluid Dynamics” - GA 681447, H2020-ERC CoG 2015 AROMA-CFD, PI G. Rozza, and INdAM-GNCS 2019–2020 projects. This work was also partially supported by the U.S. National Science Foundation through Grant No. DMS-1953535 (PI A. Quaini). A. Quaini acknowledges the support from the Radcliffe Institute for Advanced Study at Harvard University, where she has been the 2021–2022 William and Flora Hewlett Foundation Fellow.

The authors have no conflicts to disclose.

Michele Girfoglio: Investigation (equal); Methodology (equal); Validation (equal). Annalisa Quaini: Investigation (equal); Methodology (equal); Validation (equal). Gianluigi Rozza: Funding acquisition (lead); Supervision (lead).

The code created for this paper is available within open-source package Geophysical and Environmental Applications (GEA).⁷