What: A total of 47 participants from 9 countries representing 15 different oceanic numerical models met to review our current understanding of future challenges in the design of oceanic dynamical cores.
When: 17–19 September 2018
Where: Paris, France
Oceanic numerical models are used to understand and predict a wide range of processes from global paleoclimate scales to short-term prediction in estuaries and shallow coastal areas. One of the overarching challenges, and the main topic of the Community for the Numerical Modeling of the Global, Regional, and Coastal Ocean (COMMODORE) workshop is the appropriate design of the dynamical cores given the wide variety of scales of interest and their interactions with atmosphere, sea ice, biogeochemistry, and even societal processes. The construction of a dynamical core is a very long effort that takes years and decades of research and development, and requires a collaborative mixture of scientific disciplines. This work involves a significant number of fundamental choices, such as which equations to solve, which horizontal and vertical grid arrangement is adequate, which discrete algorithms allow jointly computational efficiency and sufficient accuracy, etc. Nowadays, a broad range of numerical methods are implemented in models used for realistic ocean simulations, and, owed to the advances in computational power, a meeting point has been reached between global circulation models and regional local models, such that there can be mutual benefits of a cross-fertilization between communities. This report outlines an initiative to bring together the worldwide leading researchers actively contributing to the development of oceanic model dynamical cores, such that participants could network together and focus on next challenges irrespective of target applications (regional, coastal, or global). The first COMMODORE workshop (https://commodore2018.sciencesconf.org/) was organized in Paris, France, in September 2018. In total, the participants represented 15 oceanic dynamical cores among the most widely used by the research and operational community. The motivations, topics of discussion sessions, and outcomes of the workshop are summarized below.
The ocean model developers community has had the tendency to be split depending on the target applications (global vs coastal) and on the type of horizontal grids (structured vs unstructured) and has been organized around relatively small modeling groups. However, the models have now reached such a high level of complexity that model development goes beyond the expertise of one given group and requires interactions between physicists, mathematicians, and computer scientists. In this context, this workshop aimed at gathering a community of “model oriented” researchers to foster more regular exchanges and share expertise on outstanding issues and perspectives. During this first workshop, the emphasis was on reviewing the characteristics and diversity in the formulation of oceanic models used for realistic applications as well as on outlining upcoming challenges.
EVOLUTION OF OCEANIC MODELS.
Historically, global and regional ocean models have been based on the hydrostatic primitive equations (e.g., Griffies and Adcroft 2008) discretized on structured grids using a mixture of finite-difference and finite-volume techniques for the discretization in space. The time dimension is usually treated using standard predictor–corrector or two-level approaches (e.g., Lemarié et al. 2015b). Those choices have been made because of their good compromise between simplicity, efficiency, and accuracy. In recent years, significant progress has been made for ocean modeling on unstructured grids, either via the finite volume (e.g., Chen et al. 2003; Ringler et al. 2010; Danilov et al. 2017) or finite element (e.g., Zhang et al. 2016; Korn 2017; Kärnä et al. 2018) approach. Unstructured grid models have reached an unprecedented level of maturity at least for two reasons. First, the vertical dimension is treated in a structured way compared to earlier initiatives trying to get three-dimensional unstructured meshes working. Second, a better understanding of computational modes and dispersion properties associated with a wide range of possible choices of finite-element pairs has been reached (e.g., Le Roux et al. 2007; Le Roux 2012; Eldred and Le Roux 2018). For example, Korn and Danilov (2017) have recently proposed a specific mimetic approach to control the well-known spurious mode occurring in triangular C grids (Wolfram and Fringer 2013). Unstructured grid models have been used for coastal applications for many years and they have now reached the application phase for global applications (e.g., Sidorenko et al. 2015; Petersen et al. 2018, manuscript submitted to J. Adv. Model. Earth Syst.). A long-standing concern is that computational cost per nominal grid point is generally much larger than for structured grid models, and this problem is further compounded by the absence of time refinement to locally adjust the time step to the mesh resolution when explicit time stepping is used. The test strategy presented in the next section should provide a way to quantify more rationally the difference in terms of computational costs among existing models.
Recent advances also include the development of hybrid (or generalized) vertical coordinate systems based on arbitrary Lagrangian–Eulerian (ALE) methods. The vertical distribution of Eulerian coordinate levels is predefined, whereas Lagrangian coordinate levels freely evolve with the flow. ALE methods combine the advantages of well-defined (i.e., undistorted) meshes and reduced numerical mixing, and also allow adaptation strategies (e.g., Bleck 2002; Burchard and Beckers 2004; White et al. 2009; Leclair and Madec 2011; Petersen et al. 2015). A difficulty in this case is the rezoning (also known as regridding) phase to maintain the integrity of the grid locally and globally.
A tendency in the design of the oceanic dynamical cores is the extension to the nonhydrostatic (NH) equations. The most widely used approach nowadays in oceanic models is based on the incompressible NH system solved using a pressure correction–projection method that requires the solution of a 3D Poisson equation (Lai et al. 2011; Vitousek and Fringer 2013; Voltzinger and Androsov 2016). Recently, Auclair et al. (2018) have proposed the use of the compressible nonhydrostatic ocean equations with the advantage that no global algebraic system needs to be solved to compute NH pressure anomalies, but with the disadvantage to permit acoustic modes. In this case, a specific numerical procedure is required to maintain acceptable stability of the whole code.
TEST STRATEGY AND BENCHMARK SUITE.
Given the wide variety of choices that need to be made during the development of dynamical cores and their overall complexity, it is crucial to define evaluation methods to compare the behavior of different models. Such effort has been made over the last decade by the global atmospheric community [Dynamical Core Model Intercomparison Project (DCMIP); Ullrich et al. 2012]. In particular, within DCMIP, a collection of test cases that found broad acceptance in the community has been designed and applied by a large number of modeling groups. The workshop highlighted that in the context of the oceanic community, existing test cases are scattered in the literature and not always fully documented and reproducible. The few existing examples of such effort (e.g., Ezer et al. 2002; Gerritsen et al. 2008; Ilicak et al. 2012; Soufflet et al. 2016) turned out to provide valuable feedback on the consequences of model formulations. A good test case should be easy to configure with analytical data suitable for all horizontal grids and different vertical coordinates and easy to evaluate while being relevant to test a given component of the dynamical core. The evaluation can be done either via analytical solutions (e.g., Bristeau et al. 2018), numerically converged solutions (provided that all models converge toward the same solution), or more subjectively based on an unambiguous physical understanding of the processes (e.g., Marques et al. 2017). Such a benchmark suite is also useful to motivate communication between modeling groups and also to open room for prospective approaches from applied mathematicians to highlight their effectiveness. Existing initiatives aiming at oceanic model intercomparison based on realistic simulations (e.g., Chassignet et al. 2000; Griffies et al. 2009) are generally too complex to clearly associate the observed differences to particular numerical choices. A way to evaluate numerical models in such complex configurations that could gain ground in the next few years is the uncertainty quantification (e.g., Iskandarani et al. 2016). Such an approach provides tools to characterize modeling and numerical sensitivities.
Throughout the three days of discussion, different current and future challenges have been identified. Addressing these challenges requires closer collaboration between the modeling groups.
Multiresolution strategy: Block structured versus unstructured.
On the one hand, unstructured grid models have the ability to allow variable-resolution meshes provided an efficient mesh generation tool (e.g., Engwirda 2017), for example, to adapt the resolution to follow the local Rossby radius (Hallberg 2013; Sein et al. 2017). On the other hand, structured grid models can also locally increase the resolution via nesting techniques (Debreu and Blayo 2008; Warner et al. 2010; Debreu et al. 2012) or quadtree–octree refinement (Popinet and Rickard 2007). One advantage of the nesting approach is to allow the adjustment of the time step and the physical parameters to the local resolutions while unstructured models will need scale-aware parameterizations and a specific procedure for time refinement. Other points to investigate are the impact of variable resolution on propagating waves and the optimal layout to build a multiresolution mesh. A more prospective approach could be the use of an adaptive wavelet method (Kevlahan et al. 2015).
Energy consistency and resolved/unresolved scales coupling.
Energy consistency is an important aspect for the proper interaction between resolved and parameterized scales (e.g., Burchard 2002; Bachman et al. 2017). However, as soon as a numerical core does not globally conserve energy at a discrete level (e.g., because of monotonicity enforcement, vertical remapping, or some form of upwinding), the identification of energy pathways is difficult and requires an in-depth analysis to close the energy cycle (e.g., Marsaleix et al. 2008; Eden 2016), which can be rather tedious (if not impossible) when advanced high-order numerics are used. An alternative could be to opt for an energy-conserving space and time discretization (e.g., Korn 2017; Eldred et al. 2019) with specific care (based on explicit numerical dissipation in the dynamical core and/or on parameterization of subgrid processes) to avoid instability issues of the existing approaches (Bell et al. 2017). This consistent coupling between dynamical cores and subgrid processes (known as physics–dynamics coupling) is an increasingly important topic for the building of geophysical models in general (Gross et al. 2018).
Vertical coordinates and spurious numerical mixing.
Spurious mixing (especially spurious dianeutral mixing) is a long-standing issue for oceanic dynamical cores (e.g., Griffies et al. 2000). The use of an ALE vertical coordinate is a way to mitigate this issue. Despite the fact that Gibson et al. (2017) suggest that the vertical component of spurious mixing dominates as horizontal resolution increases, it should not overshadow that many components of dynamical cores can be a source of numerical mixing (e.g., horizontal advection, time stepping, the stabilization of the mode-splitting procedure). There is still a need to better understand the implications of different choices of momentum/tracers advection schemes, rezoning, and remapping procedures on numerical mixing. Idealized test cases and efficient diagnostic tools are important to tackle this issue (e.g., Ilicak et al. 2012; Klingbeil et al. 2014). It would also be instructive to keep investigating the improvement of quasi-Eulerian vertical coordinates (e.g., Berntsen 2011). In this context schemes for the internal pressure gradient and for isoneutral tracer diffusion should also be considered (e.g., Shao et al. 2018).
Nonhydrostatic pressure contribution.
As discussed earlier, many oceanic dynamical cores have the possibility to account for nonhydrostatic effects. A difficulty is now to clearly identify under which conditions relaxing the hydrostatic assumption is necessary and which resolution is required for proper NH modeling. Another challenge is the possibility to locally account for NH effects within a primitive equations model either in the form of a superparameterization (e.g., Campin et al. 2011) or in the form of two-way nesting between coarse hydrostatic and fine nonhydrostatic meshes (e.g., Blayo and Rousseau 2016). Ultimately, further investigations of the merits and flaws of the incompressible versus compressible NH approaches would be worthwhile, also for global applications (Losch et al. 2004).
Coupling with other Earth system compartments.
Oceanic dynamical cores are often used as a component of larger coupled model systems. Coupling to several other Earth system compartments is common, such as to surface wave, sea ice, atmospheric, biogeochemical, benthic, and hydrological models. The numerical implementation of such coupling can become an issue (e.g., Lemarié et al. 2015a; Beljaars et al. 2017), especially at high coupling frequency and/or spatial resolution. More systematic analysis of the coupling stability and consistency using simplified equation sets and the design of simplified coupled test cases must be encouraged.
Vanishing layers, wetting and drying, and shock-resolving numerics.
An accurate treatment of wetting and drying is essential for coastal simulations as well as for climate simulations of under-ice-shelf cavities. At a numerical level this requires the nonnegativity of the water height and an adequate volume-conserving treatment of dry states (also known as vacuum states). Standard numerical methods used in oceanic dynamical cores do not have the ability to handle vacuum, or equivalently shocks. Instead, approaches based on some predefined minimum water depth and specific ad hoc manipulation of discrete fluxes are often used (e.g., section 5.2 in Klingbeil et al. 2018). However, there is a vast literature dedicated to the design of numerical schemes preserving positivity and able to correctly treat vacuum states that furthermore satisfy an entropy-preserving property; that is, the nonlinear solution is physically relevant even in the presence of discontinuities (e.g., Audusse et al. 2004, 2016). Considering advection–diffusion equations, for example, for tracers (temperature, salinity), there is an obvious benefit in using numerical schemes preserving the maximum principle. There could be an interest in comparing these more advanced approaches with the usual treatment adopted in dynamical cores.
The workshop gave a broad and fresh overview of existing numerical methods used in realistic ocean models as well as some examples of alternatives from the applied maths community. The participants have been enthusiastic and very positive about the possibility to sustain this type of workshop into a biennial workshop series. A collective article is currently in preparation to summarize the challenges and prospects for oceanic numerical cores across all scales. Moreover, a particular effort will be directed toward the formation of an active community of model developers with international collaborations, starting, for example, with the standardization of existing idealized test cases as the basis for model/methods intercomparison studies. The next meeting will be organized either in fall 2019 or winter 2020 in Hamburg, Germany.
The authors thank Gurvan Madec for reading and checking the content of this summary as well as all of the workshop’s participants (https://commodore2018.sciencesconf.org/resource/listeparticipants). The workshop was hosted by the French National Institute for computer science and applied mathematics (Inria), Paris, France. We gratefully acknowledge the support from Inria, from the Coastal and Regional Ocean Community Model research group (Gdr Croco), and from the French program Earth’s Fluid Envelopes (LEFE).