A dynamical core based on a discontinuous Galerkin method for higher-order finite-element sea ice modeling

The ability of numerical sea ice models to reproduce localized deformation features associated with fracture processes is key for an accurate representation of the ice dynamics and of dynamically coupled physical processes in the Arctic and Antarctic. Equally key is the capacity of these models to minimize the numerical diffusion stemming from the advection of these features to ensure that the associated strong gradients persist in time, without the need to unphysically re-inject energy for re-localization. To control diffusion and improve the approximation quality, we present a new numerical core for the dynamics of sea ice that is based on higher-order finite-element discretizations for the momentum equation and higher-order discontinuous Galerkin methods for the advection. The mathematical properties of this core are discussed, and a detailed description of an efficient shared-memory parallel implementation is given. In addition, we present different numerical tests and apply the new framework to a benchmark problem to quantify the advantages of the higher-order discretization. These tests are based on Hibler's viscous–plastic sea ice model, but the implementation of the developed framework in the context of other physical models reproducing a strong localization of the deformation is possible.

1 Introduction

Sea ice plays a critical role for the development of the Earth system, with up to 15 % of the world's oceans being covered by it at some point during the year. It contributes importantly to the global energy budget, and its high albedo keeps Arctic and Antarctic oceans cool, affecting global oceanic circulation. An accurate simulation of sea ice is therefore of importance, in particular to describe the evolution and impact of climate change. The numerical modeling of sea ice is, however, very challenging since it is characterized by nonlinear and highly localized processes.

In the present work, we develop a numerical scheme for sea ice that uses higher-order finite elements for the sea ice momentum and the advection equations and specifically aims to provide a high-fidelity discretization with small numerical diffusion and good approximation properties. We choose discontinuous Galerkin methods for the advection because they allow a Eulerian treatment of the equations of motion that is compatible with the habits of the sea ice and climate modeling community while extending naturally to high order and exhibiting limited numerical diffusion. The momentum equation will be formulated in a variational finite-element way that also allows naturally for higher-order schemes and allows for a direct coupling to the discontinuous Galerkin advection discretization. The proposed numerical scheme will form the dynamical core in the next-generation sea ice model with discontinuous Galerkin (neXtSIM-DG) that is currently under development.

Since the first extended in situ observational campaigns of the 1970s in the Arctic, sea ice has been identified as a densely fractured material in which most of the deformation is taking place locally by relative motion at the cracks with the ice between the cracks being virtually rigid (Coon et al., 1974). This relative motion of ice plates, referred to in the sea ice community as floes, translates into three main deformation processes: opening of fractures; joining along larger features called leads; and the shearing along opened fractures and the closing of leads, resulting in the formation of pressure ridges. Although highly localized around cracks, the processes play a key role in the polar ocean systems by governing the location and intensity of bio-chemical processes and the exchange of heat, mass, and momentum between the ice, ocean, and atmosphere (e.g., Marcq and Weiss, 2012; Vihma, 2014; Goosse et al., 2018; Horvat and Tziperman, 2018; Taylor et al., 2018). Importantly, the three processes also determine to a significant extent the large-scale mechanical resistance of the ice cover and hence its mobility and the overall rates of ice export out of the Arctic (e.g., Rampal et al., 2009, 2011).

Satellite remote sensing data, such as the RADARSAT Geophysical Processor System sea ice motion products which became available in the late 1990s, have allowed for the observation of these localized processes at the global scale of the Arctic Ocean. The term linear kinematic features (LKFs) was then proposed to designate the associated near-linear zones of discontinuities in the drift velocity fields. These LKFs correspond to areas with a high density of fractures in the ice cover, which strongly concentrates its deformation (Kwok, 2001). In recent years, a large number of observational analyses of sea ice deformation data, e.g., Lindsay and Stern (2003), Marsan et al. (2004), Rampal et al. (2008), Stern and Lindsay (2009), Hutchings et al. (2011), and Oikkonen et al. (2017), have fueled a race in the modeling community towards a better reproduction of LKFs in thermodynamical models, in particular with respect to their spatial and temporal statistics, e.g., Girard et al. (2011), Rampal et al. (2016), Hutter et al. (2018), Rampal et al. (2019), and Bouchat et al. (2022). Different approaches have been taken towards this goal: new mechanical (i.e., rheological) continuum models have been proposed for sea ice (Schreyer et al., 2006; Sulsky and Peterson, 2011; Girard et al., 2011; Dansereau et al., 2016; Ólason et al., 2022), the mechanical parameters of existing models have been tuned (Bouchat and Tremblay, 2017), and the spatial resolution of models has been increased (Hutter et al., 2018).

The ability to adequately reproduce LKFs in continuum sea ice models, however, raises an equally important challenge: that of keeping the very strong gradients in sea ice properties (e.g., velocity, thickness, concentration) that stem from the extreme localization of the deformation as the ice is advected by winds and ocean currents. This numerical discretization problem is, in fact, not unique to sea ice but encountered for all materials that are experiencing both highly localized deformations resulting from brittle fracturing processes and high post-fracture strains. Another important example from the geosciences is the Earth crust, where brittle processes leading to strain localization and slip coexist in faults, landslides, and volcanic edifices, e.g., Peng and Gomberg (2010) and Burov (2011). Sea ice, however, represents an extreme case as it is constantly moving and experiencing much larger relative deformations and drift velocities (about 5 and 10 cm s⁻¹ as a daily average in the winter and summer, respectively).

Several numerical approaches have been studied, and dedicated advection schemes have been developed to limit numerical diffusion in models of the Earth crust (e.g., see Zhong et al., 2015, for a review). In the sea ice modeling community, however, the treatment and, in particular, the quantification of numerical diffusion of advected gradients has received relatively little attention. Notable exceptions are the works by Flato (1993) and Huang and Savage (1998), who applied particle-in-cell methods to treat the advection of strong gradients in ice concentration and thickness not associated with LKFs but with the migration of the edge of the Arctic sea ice cover (the so-called ice edge); by Lipscomb and Hunke (2005), who used an incremental remapping to limit diffusion; by Sulsky and Peterson (2011), who introduced the material point method and tested its robustness in the context of sea ice by performing idealized convection benchmark problems; and by Danilov et al. (2015), who employed a flux-corrected Taylor–Galerkin method. neXtSIM (Rampal et al., 2016) is based on a Lagrangian model and hence completely avoids diffusion during transport, although remeshing operations are required in this framework which themselves induce some diffusion. The implementation of discontinuous Galerkin methods to treat the advection of sea ice was first proposed by Dansereau et al. (2016, 2017) and used with higher orders, with a quantification of diffusion by Dansereau et al. (2021). Mehlmann et al. (2021b) compared sea ice simulations using different meshes, mesh resolutions, and advection schemes. However, the focus of their paper was the discretization of the momentum equation, and no specific discussion of numerical diffusion was given.

Outline. The following section will introduce the basic equations and the notation used throughout the paper. We limit ourselves to the most widely used dynamical framework, which is the so-called visco-plastic rheology (Hibler, 1979), to focus on the discretization and to aid comparison to other numerical schemes in the literature. We will extend the discretization to more recently developed elasto-brittle schemes (Maxwell elasto-brittle (MEB) and brittle Bingham–Maxwell (BBM), e.g., Dansereau et al., 2016; Ólason et al., 2022) elsewhere. Section 3 details the numerical discretization of the sea ice model, including the advection and the momentum equations. Section 4 focuses on the implementation as well as on the shared-memory parallelization of the numerical model. Finally, in Sect. 5 we consider basic tests to validate the method and apply it to established benchmark problems (Mehlmann and Korn, 2021). The paper concludes with an outlook.

2 Governing equations

We denote by Ω⊂ℝ² the two-dimensional domain of the sea ice. The sea ice models we investigate consist of a momentum equation for the velocity field $\mathbf{v}:\mathrm{\Omega }\to {\mathbb{R}}^{\mathrm{2}}$ and further advection equations for tracer variables. In simple models, such as the one introduced by Hibler (1979), the tracers are usually the mean ice height $H:\mathrm{\Omega }\to [\mathrm{0},\mathrm{\infty })\subset \mathbb{R}$ and ice concentration $A:\mathrm{\Omega }\to [\mathrm{0},\mathrm{1}]\subset \mathbb{R}$. Here, we consider the following system of sea ice equations:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& {\mathit{\rho }}_{\text{ice}}H{\partial }_t\mathbf{v}=\text{div}\mathit{\sigma }+A\mathit{\tau }\left(\mathbf{v}\right)-{\mathit{\rho }}_{\mathrm{ice}}H(f_{\mathrm{c}}{\mathit{e}}_z\times \mathit{v}+g\mathrm{\nabla }{\stackrel{\mathrm{̃}}{H}}_{\mathrm{g}}),\\ & {\partial }_{t}A+\text{div}\left(\mathbf{v}A\right)=\mathrm{0},{\partial }_{t}H+\text{div}\left(\mathbf{v}H\right)=\mathrm{0}.\end{array}\end{array}$$

Here, ρ_ice is the ice density, f_ce_z×v is the Coriolis term with Coriolis parameter f_c and vertical unit vector e_z, g is the gravitational acceleration, and ${\stackrel{\mathrm{̃}}{H}}_{\mathrm{g}}$ the sea surface height. We focus on a stand-alone dynamics model without coupling to an ocean and an atmospheric model. Following Coon (1980), we approximate the surface height by the Coriolis term

$$-{\mathit{\rho }}_{\text{ice}}Hg\mathrm{\nabla }{\stackrel{\mathrm{̃}}{H}}_{\mathrm{g}}\approx {\mathit{\rho }}_{\text{ice}}H{f}_{\mathrm{c}}{\mathit{e}}_z\times {\mathbf{v}}_{\mathrm{o}},$$

where v_o is the ocean surface velocity. The forcing τ(v) is given by

$$\mathit{\tau }\left(\mathbf{v}\right)=C_{\mathrm{o}}{\mathit{\rho }}_{\mathrm{o}}\Vert {\mathbf{v}}_{\mathrm{o}}-\mathbf{v}{\Vert }_{\mathrm{2}}\cdot ({\mathbf{v}}_{\mathrm{o}}-\mathbf{v})+C_{\mathrm{a}}{\mathit{\rho }}_{\mathrm{a}}\Vert {\mathbf{v}}_{\mathrm{a}}{\Vert }_{\mathrm{2}}\cdot {\mathbf{v}}_{\mathrm{a}}.$$

The index “o” represents the ocean with the surface drag C_o, the water density ρ_o, and again the ocean surface velocity v_o, while “a” denotes the atmosphere with drag coefficient C_a, density ρ_a, and wind field v_a. We neglect turning angles and thermodynamic effects in Eq. (1). Therefore, the constraints $A\in [\mathrm{0},\mathrm{1}]$ and $H\in [\mathrm{0},\mathrm{\infty })$ are not naturally enforced by the equations but must be ensured by projections. In the following, we will use the following notation of the momentum equation (with the approximation of the surface height):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& {\mathit{\rho }}_{\text{ice}}H{\partial }_t\mathbf{v}=\text{div}\mathit{\sigma }+F\left(\mathbf{v}\right),\\ & F\left(\mathbf{v}\right)=A\mathit{\tau }\left(\mathbf{v}\right)+{\mathit{\rho }}_{\text{ice}}H{f}_{\mathrm{c}}{\mathit{e}}_z\times ({\mathbf{v}}_{\mathrm{o}}-\mathbf{v}).\end{array}\end{array}$$

Model (1) is closed by specifying a rheology, i.e., the relation between the (vertically integrated) stress σ and the strain rate ϵ,

$$\mathit{ϵ}\left(\mathbf{v}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathrm{\nabla }\mathbf{v}+\mathrm{\nabla }{\mathbf{v}}^T),{\mathit{ϵ}}^{\prime }\left(\mathbf{v}\right)=\mathit{ϵ}\left(\mathbf{v}\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\text{tr}\left(\mathit{ϵ}\right(\mathbf{v}\left)\right)I,$$

as well as the ice tracer quantities H and A (and possibly further parameters). Different rheological models have been proposed in the literature. As this paper focuses on computational questions that are largely independent of the chosen rheology, we consider the most widely used one, i.e., the viscous–plastic (VP) model proposed by Hibler (1979). It prescribes

$$\begin{array}{}\text{(3)}& \mathit{\sigma }\left(\mathbf{v}\right)=\mathrm{2}\mathit{\eta }{\mathit{ϵ}}^{\prime }\left(\mathbf{v}\right)+\mathit{\zeta }\text{div}\left(\mathbf{v}\right)I-{\displaystyle \frac{P}{\mathrm{2}}}I,\end{array}$$

with viscosities η,ζ that, using the notation introduced in Mehlmann and Richter (2017), are given by

$$\begin{array}{}\text{(4)}& \mathit{\eta }={\displaystyle \frac{\mathit{\zeta }}{e^{\mathrm{2}}}},\mathit{\zeta }={\displaystyle \frac{P_{\mathrm{0}}}{\mathrm{2}\sqrt{{\mathrm{\Delta }}_{\text{min}}^{\mathrm{2}}+\text{tr}(\mathit{ϵ}{)}^{\mathrm{2}}+\mathrm{2}{e}^{-\mathrm{2}}\cdot {\mathit{ϵ}}^{\prime }:{\mathit{ϵ}}^{\prime }}}}.\end{array}$$

Here e=2 is the eccentricity of the elliptical yield curve, Δ_min>0 is the threshold defining the transition to a viscous behavior for very small strain, P₀ is the ice strength, and P is the replacement pressure:

$$\begin{array}{}\text{(5)}& P_{\mathrm{0}}=P^{\star }\cdot H\cdot \exp (-C(\mathrm{1}-A\left)\right),P={\displaystyle \frac{\mathrm{\Delta }\left(\mathit{ϵ}\right)}{{\mathrm{\Delta }}_{\text{min}}+\mathrm{\Delta }\left(\mathit{ϵ}\right)}}\cdot P_{\mathrm{0}}.\end{array}$$

Common default values for the model parameters ${\mathit{\rho }}_{\text{ice}},{\mathit{\rho }}_{\text{a}},{\mathit{\rho }}_{\mathrm{w}},e,C,P^{\star }$ can be found in Table 1.

The VP model is highly nonlinear. Therefore, a solution with implicit methods is very challenging; see Losch et al. (2014), Mehlmann and Richter (2017), and Shih et al. (2023) for various approaches based on Newton's method. Picard iterations are also slow, and an explicit time stepping would require excessively small time steps (Ip et al., 1991). Hence, the so-called elastic–viscous–plastic (EVP) model is a widely used variant of the VP rheology (Hunke, 2001; Kimmritz et al., 2016). It adds a pseudo-elastic behavior to improve numerical performance. The constitutive law (3) is in this case given by

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{1}}{E}}\cdot {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\mathit{\sigma }+\mathit{\sigma }=\mathit{\sigma }\left(\mathbf{v}\right),\end{array}$$

where σ(v) is the VP relation given by Eq. (3). The EVP model should, however, be considered as a model different from VP since its solutions do not converge to the VP ones. An alternative variant that can be considered as a pseudo-time-stepping scheme is the mEVP scheme (Bouillon et al., 2013); see Sect. 3.4. The mEVP scheme converges to the VP solution given a sufficiently large number of iterations.

3 Higher-order finite-element discretization of the sea ice equations

In the following, we describe the discretization of the sea ice Eq. (1) in space and time using higher-order finite elements. All tracers and also the strain rate tensor ϵ and the stresses σ are discretized with a discontinuous Galerkin (dG) approach, whereas the ice velocity v is discretized using quadratic continuous finite elements.

3.1 Mesh domain

Discretizations of the sea ice equations are typically used within a coupled Earth system model. One consequence is that the time step of the numerical sea ice model is not only determined by the desired accuracy and stability considerations but also constrained by the atmospheric and oceanic components of the Earth system model.

By Δt we denote the time step of the sea ice equations. Although dynamic time discretizations with varying step sizes are possible, we will only consider uniform time steps with Δt_n=Δt for all steps n. The time mesh is hence given by

$$\begin{array}{}\text{(7)}& t_{\mathrm{0}}<t_{\mathrm{1}}<t_{\mathrm{2}}<\mathrm{\dots }<t_N=T,\mathrm{\Delta }t:=t_n-t_{n-\mathrm{1}}.\end{array}$$

Assuming, for example, the time step size Δt∈240 s and ice velocities $|\mathbf{v}{|}_{\mathrm{\infty }}\le \mathrm{1}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$, explicit time stepping will be stable for mesh sizes up to a resolution of

$$\begin{array}{}\text{(8)}& \mathrm{\Delta }x\approx C_r\cdot \Vert \mathbf{v}\Vert \cdot \mathrm{\Delta }t\approx C_r\cdot \mathrm{250}\mathrm{m},\end{array}$$

where C_r is a constant that depends on the degree r of the time-stepping scheme. The factor C_r scales like $C_r\approx \mathrm{2}r+\mathrm{1}$ (Chalmers and Krivodonova, 2020). Hence, for a dG(2) method and time-stepping scheme of balanced order with r=2, the minimum mesh element size should be larger than 2 km if a time step of Δt=240 s is used. Usually, limitations due to the CFL condition are more relevant in the ocean model, where velocities are higher. For our higher-order sea ice model the CFL condition might become important due to its scaling with the polynomial degree.

Figure 1 Parametric mesh. Each element $T_{i,j}\in {\mathcal{T}}_h$ (on the right) arises from the a mapping ${\mathbf{T}}_{i,j}:\widehat{T}\to T_{i,j}$ from the reference element $\widehat{T}=(\mathrm{0},\mathrm{1}{)}^{\mathrm{2}}$ (on the left). The mesh elements T_i,j are general quadrilaterals such that the mappings T are bilinear polynomials. The edges $e_{i,j}^x$ and $e_{i,j}^y$ are straight lines.

For the spatial discretization, we employ a parametric finite-element mesh of the domain Ω. We base the discretization on quadrilateral meshes 𝒯_h (as opposed to triangular ones). The meshes are topologically fully regular but geometrically distorted and consist of nodes x_i,j for $i=\mathrm{0},\mathrm{\dots },N_x$, $j=\mathrm{0},\mathrm{\dots },N_y$ and elements $T_{i,j}=({\mathit{x}}_{i-\mathrm{1},j-\mathrm{1}},{\mathit{x}}_{i,j-\mathrm{1}},{\mathit{x}}_{i-\mathrm{1},j},{\mathit{x}}_{i,j})$ for $i=\mathrm{1},\mathrm{\dots },N_x$, $j=\mathrm{1},\mathrm{\dots },N_y$, horizontal edges $e_{i,j}^x=({\mathit{x}}_{i-\mathrm{1},j},{\mathit{x}}_{i,j})$ for $i=\mathrm{1},\mathrm{\dots },N_x$, $j=\mathrm{0},\mathrm{\dots },N_y$ and vertical edges $e_{i,j}^i=({\mathit{x}}_{i,j-\mathrm{1}},{\mathit{x}}_{i,j})$,

$$\begin{array}{}\text{(9)}& {\mathcal{T}}_h=\mathit{\{}{x}_{i,j},T_{i,j},e_{i,j}^x,e_{i,j}^y\mathit{\}},\end{array}$$

where $N_x,N_y\in \mathbb{N}$ denote the number of elements in the x and y direction. See Fig. 1 for an illustration. The nodes are lexicographically ordered; i.e., $k=i+(N_x+\mathrm{1})j$ is the consecutive index. Each geometric mesh element T_i,j can be defined via a mapping

$${\mathbf{T}}_{i,j}:\widehat{T}:=(\mathrm{0},\mathrm{1}{)}^{\mathrm{2}}\mapsto T_{i,j}$$

from a unique reference element $\widehat{T}$ using the bilinear polynomial

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\mathbf{T}}_{i,j}\left(\widehat{\mathbf{x}}\right)& =(\mathrm{1}-{\widehat{\mathbf{x}}}_{\mathrm{1}})(\mathrm{1}-{\widehat{\mathbf{x}}}_{\mathrm{2}}){\mathbf{x}}_{i-\mathrm{1},j-\mathrm{1}}+{\widehat{\mathbf{x}}}_{\mathrm{1}}(\mathrm{1}-{\widehat{\mathbf{x}}}_{\mathrm{2}}){\mathbf{x}}_{i,j-\mathrm{1}}\\ & +(\mathrm{1}-{\widehat{\mathbf{x}}}_{\mathrm{1}}){\widehat{\mathbf{x}}}_{\mathrm{2}}{\mathbf{x}}_{i-\mathrm{1},j}+{\widehat{\mathbf{x}}}_{\mathrm{1}}{\widehat{\mathbf{x}}}_{\mathrm{2}}{\mathbf{x}}_{i,j};\end{array}\end{array}$$

see again Fig. 1. On each edge $e_{i,j}^x$ and $e_{i,j}^y$, we consider one unit normal vector. Its orientation arises from mapping the unit normal vectors ${\widehat{e}}^x=(\mathrm{1},\mathrm{0}{)}^T$ and ${\widehat{e}}^y=(\mathrm{0},\mathrm{1}{)}^T$ of the reference element to the edges of 𝒯_h.

3.2 Finite-element spaces and degrees of freedom

We use continuous (cG) and discontinuous (dG) finite elements for the discretization of the momentum equation as well as the constitutive equations and advection problems, respectively. On the reference element $\widehat{T}$, we define two sets of basis functions. The dG basis functions that we employ are given by

$$\begin{array}{}\text{(11)}& \begin{array}{l}{\mathrm{\Psi }}_{\mathrm{1}}\left(\widehat{\mathbf{x}}\right):=\mathrm{1},\\ {\mathrm{\Psi }}_{\mathrm{2}}\left(\widehat{\mathbf{x}}\right):={\widehat{\mathbf{x}}}_{\mathrm{1}}-\frac{\mathrm{1}}{\mathrm{2}},\\ {\mathrm{\Psi }}_{\mathrm{3}}\left(\widehat{\mathbf{x}}\right):={\widehat{\mathbf{x}}}_{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}},\\ {\mathrm{\Psi }}_{\mathrm{4}}\left(\widehat{\mathbf{x}}\right):=\left({\widehat{\mathbf{x}}}_{\mathrm{1}}-\frac{\mathrm{1}}{\mathrm{2}}\right)\left({\widehat{\mathbf{x}}}_{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}}\right),\\ {\mathrm{\Psi }}_{\mathrm{5}}\left(\widehat{\mathbf{x}}\right):={\left({\widehat{\mathbf{x}}}_{\mathrm{1}}-\frac{\mathrm{1}}{\mathrm{2}}\right)}^{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{12}},\\ {\mathrm{\Psi }}_{\mathrm{6}}\left(\widehat{\mathbf{x}}\right):={\left({\widehat{\mathbf{x}}}_{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}}\right)}^{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{12}},\\ {\mathrm{\Psi }}_{\mathrm{7}}\left(\widehat{\mathbf{x}}\right):=\left({\widehat{\mathbf{x}}}_{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}}\right)\left({\left({\widehat{\mathbf{x}}}_{\mathrm{1}}-\frac{\mathrm{1}}{\mathrm{2}}\right)}^{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{12}}\right),\\ {\mathrm{\Psi }}_{\mathrm{8}}\left(\widehat{\mathbf{x}}\right):=\left({\widehat{\mathbf{x}}}_{\mathrm{1}}-\frac{\mathrm{1}}{\mathrm{2}}\right)\left({\left({\widehat{\mathbf{x}}}_{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}}\right)}^{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{12}}\right).\end{array}\end{array}$$

These basis functions are orthogonal, i.e., ${\int }_{\widehat{T}}{\mathrm{\Psi }}_i{\mathrm{\Psi }}_j\text{d}\mathbf{x}={\mathit{\delta }}_{ij}$. Second, we use the degree r tensor product Lagrange finite-element basis functions

$$\begin{array}{}\text{(12)}& {\mathrm{\Phi }}_{(r+\mathrm{1})l+k}^{\left(r\right)}\left(\widehat{\mathbf{x}}\right):={\mathit{\xi }}_k^{\left(\mathrm{1}\right)}\left({\mathbf{x}}_{\mathrm{1}}\right){\mathit{\xi }}_{\mathrm{l}}^{\left(\mathrm{1}\right)}\left({\mathbf{x}}_{\mathrm{2}}\right),k,l=\mathrm{1},\mathrm{\dots },r+\mathrm{1}.\end{array}$$

The one-dimensional basis functions (for r=1 and r=2) are given by

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\mathit{\xi }}_{\mathrm{1}}^{\left(\mathrm{1}\right)}\left(\widehat{x}\right)& :=\mathrm{1}-\widehat{x},{\mathit{\xi }}_{\mathrm{2}}^{\left(\mathrm{1}\right)}\left(\widehat{x}\right):=\widehat{x},\\ {\mathit{\xi }}_{\mathrm{1}}^{\left(\mathrm{2}\right)}\left(\widehat{x}\right)& :=(\mathrm{1}-\widehat{x})(\mathrm{1}-\mathrm{2}\widehat{x}),{\mathit{\xi }}_{\mathrm{2}}^{\left(\mathrm{2}\right)}\left(\widehat{x}\right):=\mathrm{4}\widehat{x}(\mathrm{1}-\widehat{x}),\\ {\mathit{\xi }}_{\mathrm{3}}^{\left(\mathrm{2}\right)}\left(\widehat{x}\right)& :=\widehat{x}(\mathrm{2}\widehat{x}-\mathrm{1}).\end{array}\end{array}$$

All basis functions are mapped from the reference element $\widehat{T}$ onto the mesh elements of 𝒯_h.

We define continuous finite-element spaces $V_h^{\left(r\right)}$, where $r=\mathrm{1},\mathrm{2}$ is the degree, and spaces $W_h^{\left(s\right)}$ are associated with the discontinuous finite elements, where $s=\mathrm{1},\mathrm{2},\mathrm{\dots }$ is the number of local basis functions,

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& V_h^{\left(r\right)}=\mathit{\{}\mathit{\varphi }\in C(\overline{\mathrm{\Omega }}):\mathit{\varphi }{|}_T={\mathrm{\Phi }}^{\left(r\right)}\circ {\mathbf{T}}_T^{-\mathrm{1}}\mathit{\}}\\ & W_h^{\left(s\right)}=\mathit{\{}\mathit{\psi }\in L^{\mathrm{2}}(\mathrm{\Omega }):\mathit{\psi }{|}_T={\mathrm{\Psi }}^{\circ }{\mathbf{T}}_T^{-\mathrm{1}}\mathit{\}},\end{array}\end{array}$$

where Ψ and Φ^(r) are linear combinations of reference basis functions according to Eqs. (11) and (12), respectively. Locally on each mesh element T∈𝒯_h, the tracer $H_h\in W_h^{\left(s\right)}$ and velocity ${\mathbf{v}}_h\in V_h^{\left(r\right)}$ are therefore described by the linear combinations of the basis functions:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{\widehat{H}}_T\left(\widehat{\mathbf{x}}\right)& :=H_h\left({\mathbf{T}}_T\right(\widehat{\mathbf{x}}\left)\right)=\sum _{j=\mathrm{1}}^s{H}_{T,j}{\mathrm{\Psi }}_j\left(\widehat{\mathbf{x}}\right),{\widehat{\mathbf{v}}}_T\left(\widehat{\mathbf{x}}\right)\\ & :={\mathbf{v}}_h\left({\mathbf{T}}_T\right(\widehat{\mathbf{x}}\left)\right)=\sum _{j=\mathrm{1}}^{N_{\mathrm{loc}}^{\text{cG}}}{\mathbf{v}}_{T,j}{\mathrm{\Phi }}_k\left(\widehat{\mathbf{x}}\right),\end{array}\end{array}$$

where $N_{\mathrm{loc}}^{\text{cG}}=(r+\mathrm{1}{)}^{\mathrm{2}}$ is the local number of unknowns in each element. An analogous representation holds for the second tracer. Finally, by ${\left(\cdot ,\cdot \right)}_T$ and ${〈\cdot ,\cdot 〉}_e$ we denote L² scalar products,

$${\left(\mathit{\varphi },\mathit{\psi }\right)}_T=\underset{T}{\int }\mathit{\varphi }\left(\mathbf{x}\right)\mathit{\psi }\left(\mathbf{x}\right)\text{d}\mathbf{x},{〈\mathit{\varphi },\mathit{\psi }〉}_e=\underset{e}{\int }\mathit{\varphi }\left(\mathbf{x}\right)\mathit{\psi }\left(\mathbf{x}\right)\text{d}s,$$

on the elements T_i,j and the edges $e_{i,j}^{\mathit{\{}x,y\mathit{\}}}$, respectively.

3.3 Discontinuous Galerkin advection scheme

We begin by describing the discretization of the advection equation

$${\partial }_{t}H+\text{div}\left(\mathbf{v}H\right)=\mathrm{0}$$

for the tracer $H:\mathrm{\Omega }\to \mathbb{R}$. We follow the notation of Di Pietro and Ern (2012, Chap. 3).

The temporal discretization will be by explicit Runge–Kutta schemes of order 1, 2, or 3, and in space we choose $H_h\in W_h^{\left(s\right)}$; see Eqs. (14) and (15). The discretization is based on the standard upwind formulation:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}\sum _{T\in {\mathcal{T}}_h}{\partial }_t{\left(H_h,\mathit{\psi }\right)}_T& -{\left(H_h\mathbf{v},\mathrm{\nabla }\mathit{\psi }\right)}_T\\ & +\sum _{e\in {\mathcal{T}}_h}{〈\mathit{\left\{}\mathit{\right\{}{H}_h\mathit{\left\}}\mathit{\right\}},\mathbf{v}\cdot {\mathit{n}}_e\left[\right[\mathit{\psi }\left]\right]〉}_e\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{〈|\mathbf{v}\cdot {\mathit{n}}_e|\cdot \left[\right[H_h\left]\right],\left[\right[\mathit{\psi }\left]\right]〉}_e=\mathrm{0}.\end{array}\end{array}$$

By {{H_h}}|_e we denote the average of the dG function H_h on an edge $e=\partial T_{\mathrm{1}}\cap \partial T_{\mathrm{2}}$ between the two elements T₁,T₂ and by [[H_h]]|_e the jump over this edge, i.e.,

$$\mathit{\left\{}\mathit{\right\{}{H}_h\mathit{\left\}}\mathit{\right\}}{|}_e={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(H_h{|}_{T_{\mathrm{1}}}+H_h{|}_{T_{\mathrm{2}}}),\left[\right[H_h\left]\right]{|}_e=H_h{|}_{T_{\mathrm{1}}}-H_h{|}_{T_{\mathrm{2}}}.$$

The upwind scheme can be written in matrix-vector notation as

$$\mathbf{M}{\partial }_t{H}_h=\mathbf{A}\left({\mathbf{v}}_h\right)H_h,$$

where M is the dG mass matrix in $W_h^{\left(s\right)}$, which is block diagonal with blocks of size s×s, and where A(v_h) gathers all remaining terms of Eq. (16) which are all linear in H_h. The equation is discretized in time by standard explicit Runge–Kutta methods on the advection time mesh in Eq. (7).

For dG(0) with space $W_h^{\left(\mathrm{1}\right)}$, the discretization is equivalent to the usual finite-volume upwind scheme since the per element term (H_hv,∇ψ) vanishes for all $\mathit{\psi }\in W_h^{\left(\mathrm{1}\right)}$ as ψ|_T is constant on T. The advantage of using higher-order methods will become clear in Sect. 5.3.3.

3.4 Discretizing the momentum equation

The coupled advection and momentum equation system in Eq. (1) is decoupled in a partitioned iteration by performing the advection step and then solving the momentum equation. The momentum equation is approximated with an mEVP solver, which can be considered as a pseudo time-stepping scheme for the implicit backward Euler discretization of the VP formulation (e.g., see Lemieux et al., 2012; Bouillon et al., 2013). We introduce the iterates ${\mathbf{v}}_n^{\left(p\right)}$ and ${\mathit{\sigma }}_n^{\left(p\right)}$ for $p=\mathrm{0},\mathrm{1},\mathrm{\dots },N_{\text{mEVP}}$ with ${\mathbf{v}}_n^{\left(\mathrm{0}\right)}:={\mathbf{v}}_{n-\mathrm{1}}$ and ${\mathit{\sigma }}_n^{\left(\mathrm{0}\right)}:={\mathit{\sigma }}_{n-\mathrm{1}}$, in which case the update can be written as

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& (\mathrm{1}+\mathit{\alpha }){\mathit{\sigma }}_n^{\left(p\right)}=\mathit{\alpha }{\mathit{\sigma }}_n^{(p-\mathrm{1})}+\mathit{\sigma }\left({\mathbf{v}}_n^{(p-\mathrm{1})}\right),\\ & \left(\right(\mathrm{1}+\mathit{\beta }){\mathit{\rho }}_{\mathrm{ice}}{H}_n+\mathrm{\Delta }t{A}_n{C}_{\mathrm{o}}{\mathit{\rho }}_{\mathrm{o}}\Vert {\mathbf{v}}_{\mathrm{o}}-{\mathbf{v}}_n^{(p-\mathrm{1})}{\Vert }_{\mathrm{2}}){\mathbf{v}}_n^{\left(p\right)}\\ & ={\mathit{\rho }}_{\mathrm{ice}}{H}_n({\mathbf{v}}_{n-\mathrm{1}}+\mathit{\beta }{\mathbf{v}}_n^{(p-\mathrm{1})})\\ & +\mathrm{\Delta }t(\text{div}{\mathit{\sigma }}_n^{\left(p\right)}+\stackrel{\mathrm{̃}}{F}({\mathbf{v}}_n^{(p-\mathrm{1})}\left)\right).\end{array}\end{array}$$

The forcing term $F\left({\mathbf{v}}_n^{\left(p\right)}\right)$ in Eq. (2) is split into explicit and implicit parts. The ocean forcing term is considered implicitly, which helps the stability of the scheme, and the remaining explicit terms on the right-hand side are

$$\begin{array}{}\text{(18)}& \begin{array}{rl}\stackrel{\mathrm{̃}}{F}\left({\mathbf{v}}_n^{(p-\mathrm{1})}\right)& :=A_n(C_{\mathrm{o}}{\mathit{\rho }}_{\mathrm{o}}\Vert {\mathbf{v}}_{\mathrm{o}}-{\mathbf{v}}_n^{(p-\mathrm{1})}{\Vert }_{\mathrm{2}}\cdot {\mathbf{v}}_{\mathrm{o}}+C_{\mathrm{a}}{\mathit{\rho }}_{\mathrm{a}}\Vert {\mathbf{v}}_{\mathrm{a}}{\Vert }_{\mathrm{2}}\cdot {\mathbf{v}}_{\mathrm{a}})\\ & +{\mathit{\rho }}_{\text{ice}}{H}_n{f}_{\mathrm{c}}{\mathit{e}}_z\times ({\mathbf{v}}_{\mathrm{o}}-{\mathbf{v}}_n^{(p-\mathrm{1})}).\end{array}\end{array}$$

The parameters α and β in Eq. (17) control the stability but also the speed of convergence of the mEVP iteration to the VP solution, whereas the number of steps N_mEVP controls the accuracy. We refer the reader to Kimmritz et al. (2016) for a discussion on this.

A mixed finite-element approach is used for the spatial discretization of Eqs. (17)–(18) with continuous finite elements for the momentum equation and discontinuous ones for the stress update. This yields

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& (\mathrm{1}+\mathit{\alpha })({\mathit{\sigma }}_n^{\left(p\right)},{\mathrm{\Psi }}_h)=\mathit{\alpha }({\mathit{\sigma }}_n^{(p-\mathrm{1})},{\mathrm{\Psi }}_h)+\left(\mathit{\sigma }\right({\mathbf{v}}_n^{(p-\mathrm{1})}),{\mathrm{\Psi }}_h),\\ & \left(\right((\mathrm{1}+\mathit{\beta }){\mathit{\rho }}_{\mathrm{ice}}{H}_n+\mathrm{\Delta }t{A}_n{C}_{\mathrm{o}}{\mathit{\rho }}_{\mathrm{o}}\Vert {\mathbf{v}}_{\mathrm{o}}-{\mathbf{v}}_n^{(p-\mathrm{1})}{\Vert }_{\mathrm{2}}){\mathbf{v}}_n^{\left(p\right)},{\mathrm{\Phi }}_h)\\ & =\left({\mathit{\rho }}_{\mathrm{ice}}{H}_n\right({\mathbf{v}}_{n-\mathrm{1}}+\mathit{\beta }{\mathbf{v}}_n^{(p-\mathrm{1})})+\mathrm{\Delta }t\stackrel{\mathrm{̃}}{F}({\mathbf{v}}_n^{(p-\mathrm{1})}),{\mathrm{\Phi }}_h)\\ & -\mathrm{\Delta }t({\mathit{\sigma }}_n^{\left(p\right)},\mathrm{\nabla }{\mathrm{\Phi }}_h)\end{array}\end{array}$$

for test functions

$$\begin{array}{}\text{(20)}& \begin{array}{rl}\mathrm{\Psi }\in {\mathbf{W}}_h& :=[W_h^{\left(s\right)}{]}^{\mathrm{2}\times \mathrm{2},sym}:=\mathit{\{}\mathit{\sigma }\in L^{\mathrm{2}}(\mathrm{\Omega }{)}^{\mathrm{2}\times \mathrm{2}},\\ & \mathit{\sigma }={\mathit{\sigma }}^T,{\mathit{\sigma }}_{ij}\in W_h^{\left(s\right)},i,j=\mathrm{1},\mathrm{2}\mathit{\}},\\ {\mathrm{\Phi }}_h\in {\mathbf{V}}_h& :=[V_h^{\left(r\right)}{]}^{\mathrm{2}}.\end{array}\end{array}$$

Compatibility of the velocity and stress spaces is important for the stability of the coupled iteration; see Sect. 5.2.3 for an example of possible instabilities. Stress spaces that are too small do not provide sufficient control of the term (σ(v),Ψ) in Eq. (19) and lead to oscillatory stresses; see the upper right plot in Fig. 13. The problem is related to the control of the energy, and in a mixed formulation the spaces V_h and W_h must in particular satisfy the Babuška–Brezzi condition; see, e.g., Ern and Guermond (2021, see, e.g., Theorem 49.13) for well posedness. For a simplified linear equation, this condition would mean that there exists a constant γ>0 such that

$$\begin{array}{}\text{(21)}& \underset{\mathbf{\Phi }\in {\mathbf{V}}_h}{inf}\underset{\mathbf{\Psi }\in {\mathbf{W}}_h}{sup}{\displaystyle \frac{(\mathbf{\Psi },\mathrm{\nabla }\mathbf{\Phi }{)}_{\mathrm{\Omega }}}{\Vert \mathrm{\nabla }\mathbf{\Phi }\Vert \cdot \Vert \mathbf{\Psi }\Vert }}\ge \mathit{\gamma }>\mathrm{0}.\end{array}$$

This condition can easily be satisfied if for every v_h∈V_h from the cG velocity space it holds that

$$\begin{array}{}\text{(22)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathrm{\nabla }{\mathbf{v}}_h+\mathrm{\nabla }{\mathbf{v}}_h^T)\in {\mathbf{W}}_h=[W_h^{\left(s\right)}{]}^{\mathrm{2}\times \mathrm{2},\mathrm{sym}}.\end{array}$$

Then, for any Φ=v_h we choose Ψ as Eq. (22) and get, using the symmetry of the inner product,

$$\begin{array}{rl}{\displaystyle \frac{(\mathrm{\nabla }{\mathbf{v}}_h+\mathrm{\nabla }{\mathbf{v}}_h^T,\mathrm{\nabla }{\mathbf{v}}_h{)}_{\mathrm{\Omega }}}{\Vert \mathrm{\nabla }{\mathbf{v}}_h+\mathrm{\nabla }{\mathbf{v}}_h^T{\Vert }_{\mathrm{\Omega }}\Vert \mathrm{\nabla }{\mathbf{v}}_h{\Vert }_{\mathrm{\Omega }}}}& ={\displaystyle \frac{(\mathrm{\nabla }{\mathbf{v}}_h+\mathrm{\nabla }{\mathbf{v}}_h^T,\frac{\mathrm{1}}{\mathrm{2}}(\mathrm{\nabla }{\mathbf{v}}_h+\mathrm{\nabla }{\mathbf{v}}_h^T){)}_{\mathrm{\Omega }}}{\Vert \mathrm{\nabla }{\mathbf{v}}_h+\mathrm{\nabla }{\mathbf{v}}_h^T{\Vert }_{\mathrm{\Omega }}\Vert \mathrm{\nabla }{\mathbf{v}}_h{\Vert }_{\mathrm{\Omega }}}}\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\Vert \mathrm{\nabla }{\mathbf{v}}_h+\mathrm{\nabla }{\mathbf{v}}_h^T{\Vert }_{\mathrm{\Omega }}}{\Vert \mathrm{\nabla }{\mathbf{v}}_h{\Vert }_{\mathrm{\Omega }}}}\ge {\displaystyle \frac{c_K}{\mathrm{2}}},\end{array}$$

where c_K>0 is the constant of Korn's inequality (Ern and Guermond, 2021, Theorem 42.9 and 42.10). We therefore require that the spaces W_h and V_h always allow for choosing the stress test function Ψ∈W_h as the symmetric velocity gradient, Eq. (22). To be precise, the degree s has to be chosen such that the symmetric gradient of the discrete velocity is part of the stress space. On quadrilateral elements, the continuous finite-element basis is not the pure polynomial basis P^(r), but it includes the additional mixed terms xy for r=1 and $x^{\mathrm{2}}y,x{y}^{\mathrm{2}},x^{\mathrm{2}}{y}^{\mathrm{2}}$ for r=2. Hence, the gradient space must also be enriched. For linear elements with r=1 the condition in Eq. (22) requires s=3, and for quadratic velocities with r=2 we must take s=8 in Eqs. (19)–(20). This update involves the inversion of the mass matrix of $W_h^{\left(s\right)}$. The matrix is block diagonal with block size s×s so that in the cG(2) case with s=8 the costs for the inversion are substantial. Section 4.2 describes our approach for an efficient implementation.

The momentum equation is discretized with continuous finite elements in the discrete space V_h. All zero-order terms in the momentum equation, Eq. (19), are evaluated node-wise and no integration is required. Adding the stress, however, requires integration and inversion of the mass matrix of V_h. To avoid the inversion, we use mass lumping. The evaluation of the momentum equation's right-hand side in Eq. (19) then becomes

$$\begin{array}{}\text{(23)}& \begin{array}{rl}{\mathbf{v}}_{n,i}^{\left(p\right)}& =(\mathrm{1}+\mathit{\beta }){\mathit{\rho }}_{\mathrm{ice}}{H}_{n,i}+\mathrm{\Delta }t{A}_{n,i}{C}_{\mathrm{o}}{\mathit{\rho }}_{\mathrm{o}}\Vert {\mathbf{v}}_{o,i}-{\mathbf{v}}_{n,i}^{(p-\mathrm{1})}{\Vert }_{\mathrm{2}}{)}^{-\mathrm{1}}\\ & \times \left({\mathit{\rho }}_{\mathrm{ice}}{H}_{n,i}\right({\mathbf{v}}_{n-\mathrm{1},i}+\mathit{\beta }{\mathbf{v}}_{n,i}^{(p-\mathrm{1})})+\mathrm{\Delta }t\stackrel{\mathrm{̃}}{F}({\mathbf{v}}_{n,i}^{(p-\mathrm{1})})\\ & -{\mathbf{M}}_{l,ii}^{-\mathrm{1}}\mathrm{\Delta }t({\mathit{\sigma }}_n^{\left(p\right)},\mathrm{\nabla }{\mathrm{\Phi }}_i{)}_{\mathrm{\Omega }}),\\ & i=\mathrm{1},\mathrm{\dots },N_{\mathrm{cG}},\end{array}\end{array}$$

where M_l is the lumped mass matrix in the cG velocity space. The implicit terms are handled analogously. The integration of the stresses against the gradient of the test function is a non-local operation coupling adjacent degrees of freedom. All other operations, like computing $\stackrel{\mathrm{̃}}{F}\left({\mathbf{v}}_n^{(p-\mathrm{1})}\right)$, are fully decoupled and can be processed node-wise in parallel.

Remark 1 (optimality of velocity–stress discretization).

On triangular meshes, $V_h^{\left(r\right)}$ would be the spaces of piecewise polynomials of degree r. In this case, the optimal dG stress space would be the space of piecewise constants, i.e., $W_h^{\left(\mathrm{1}\right)}$ in our notation in the case r=1 and the space of piecewise linear $W_h^{\left(\mathrm{3}\right)}$ for r=2. The choice given in Eqs. (19)–(20) appears highly inefficient as the local basis has 3 instead of 1 unknowns for r=1 and 8 instead of 3 for r=2. However, a triangular mesh with the same number of velocity unknowns has twice the number of elements as a quadrilateral mesh. Hence, r=1 has 2 unknowns on triangles and r=2 brings 6 unknowns compared to s=3 and s=8 in the case of quadrilateral meshes. This means that the difference in effort is less dramatic than it appears on first sight.

Figure 2 The numerical quadrature nodes ${\widehat{\mathit{\chi }}}_q$ are defined on the reference element $\widehat{T}=(\mathrm{0},\mathrm{1}{)}^{\mathrm{2}}$ and mapped to the real mesh elements T∈𝒯_h via ${\mathit{\chi }}_q:=\mathbf{T}\left({\widehat{\mathit{\chi }}}_q\right)$. We show a 2-point Gauss rule (of degree 4) with 2 points on each edge and 4 points in the element.

3.5 Numerical quadrature

In the parametric finite-element setup, all integrals appearing in the advection scheme in Eq. (16) and the weak formulation of the mEVP iteration in Eq. (17) must be evaluated on the reference element $\widehat{T}$ and, in case of the upwind scheme, also on the reference edge $\widehat{e}=(\mathrm{0},\mathrm{1})$ since the basis functions are defined on $\widehat{T}$. For the different terms of Eq. (16) it holds

$$\begin{array}{}\text{(24)}& \begin{array}{rl}& {\partial }_t(H_h,\mathit{\psi }{)}_T=(J_T{\partial }_t{\widehat{H}}_h,\mathrm{\Psi }{)}_{\widehat{T}}\\ & (H_h\mathbf{v},\mathrm{\nabla }\mathit{\psi }{)}_T=({\widehat{H}}_h\widehat{\mathbf{v}},J_T\widehat{\mathrm{\nabla }}{\mathbf{T}}_T^{-T}\widehat{\mathrm{\nabla }}\mathrm{\Psi }{)}_{\widehat{T}}\\ \langle & \mathit{\left\{}\mathit{\right\{}{H}_h\mathit{\left\}}\mathit{\right\}},(\mathbf{v}\cdot {\mathit{n}}_e)\left[\right[\mathit{\psi }\left]\right]{\rangle }_e=\langle J_e\mathit{\left\{}\mathit{\right\{}{\widehat{H}}_h\mathit{\left\}}\mathit{\right\}},\widehat{(\mathbf{v}\cdot {\mathit{n}}_e)}\left[\right[\mathrm{\Psi }\left]\right]{\rangle }_{\widehat{e}}\\ & {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\langle |\mathbf{v}\cdot {\mathit{n}}_e|\cdot \left[\right[H_h\left]\right],\left[\right[\mathit{\psi }\left]\right]{\rangle }_e={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\langle J_e\left|\widehat{\mathbf{v}\cdot {\mathit{n}}_e}\right|\cdot \left[\right[{\widehat{H}}_h\left]\right],\left[\right[\mathrm{\Psi }\left]\right]{\rangle }_{\widehat{e}},\end{array}\end{array}$$

where ${\widehat{H}}_h$ and $\widehat{\mathbf{v}}$ are the functions on the reference element that by Eq. (15) are associated with H_h and v on the element T and analogously for the edge terms. The reference element map T_T-dependent terms in Eq. (24) are the Jacobian $\widehat{\mathrm{\nabla }}{\mathbf{T}}_T:\widehat{T}\to {\mathbb{R}}^{\mathrm{2}\times \mathrm{2}}$ and its determinant $J_T=\text{det}\left(\widehat{\mathrm{\nabla }}{\widehat{\mathbf{T}}}_T\right)$. Since T_T is bilinear (and not linear), the Jacobian and its determinant are not element-wise constant. However, on the reference edges $\widehat{e}$, T_T is linear such that e∈𝒯_h values are straight and hence $J_e=\left|e\right|$ as the reference element has an edge length of 1.

The integrals in Eq. (24) are approximated by Gaussian quadrature. For dG(r) ($r=\mathrm{0},\mathrm{1},\mathrm{2}$) we use r+1 quadrature points on the edge and (r+1)² quadrature points within the elements; see Fig. 2 for an example with 2 points on the edges and 2×2 points within the element.

Implementation details are described in Sect. 3.1. Evaluation of the terms in Eq. (24) is numerically costly, mostly due to the evaluation of the map T_T, the Jacobian $\widehat{\mathrm{\nabla }}{\mathbf{T}}_T$, its inverse, and the determinant of the Jacobian.

4 Efficient parallelizable implementation

In the following paragraphs, we will describe the C++ implementation of the higher-order discretization. A hybrid parallelization approach consisting of distributed memory MPI splitting and local shared-memory OpenMP realization is considered. The data are structured such that the implementation also allows us to run modules on a GPU.

MPI parallelization builds on a domain decomposition that splits the complete mesh into a balanced number of rectangular subdomains such that the average number of ice-covered elements for each domain is comparable. Each parallel task then operates on a subdomain that is topologically structured into $N_{\mathrm{el}}:=N_x\times N_y$ elements such as described in Sect. 3.1.

4.1 Implementation of continuous and discontinuous finite elements

We start by describing the handling of the data, i.e., the cG and dG vectors for each MPI task that is responsible for one topologically rectangular mesh 𝒯_h consisting of N_x×N_y elements. A dG vector $A_h\in W_h^{\left(s\right)}$ has s unknowns on each of the $N_{\mathrm{el}}=N_x\cdot N_y$ elements, and we store such a vector as a $\mathbf{A}\in {\mathbb{R}}^{N_{\mathrm{el}}\times s}$ matrix. The implementation is based on Eigen (Guennebaud et al., 2010), a C++ library for linear algebra that heavily relies on C++ templates. In the code, the vector is represented as

Matrix<FloatType, Dynamic, s,
           RowMajor> DGVector<s> A;

The first dimension of A (number of elements) is dynamic and determined at runtime, which allows us to flexibly handle different subdomain sizes. The second dimension, i.e., the number of components, has degree s and is determined at compile time. This allows for vectorized SIMD processing of computations; see Sect. 5.3.2 for a numerical demonstration.

To provide one example of a frequently used operation, we explain the restriction of a dG(1) function $A_h\in W_h^{\left(\mathrm{3}\right)}$ (with three local unknowns on T) from an element T to one of its edges $e\in \partial T$. Let T have the element index $i\in \mathit{\{}\mathrm{1},\mathrm{\dots },N_{\mathrm{el}}\mathit{\}}$ and let $e=e_i^x$ be the lower edge in the notation of Eq. (9). Then the restriction to the lower edge is realized as

Vector<FloatType, 2> l_edge
(const DGVector<3>& A, size_t i)
{ return Vector<FloatType, 2> a_e(
{A(i,0)-0.5 * A(i,2), A(i,1)}); }

Since the restriction does not depend on the specific element T∈𝒯_h, the relations are implemented for the four edges and the different choices of dG spaces, i.e., for the number of local basis functions, using template specializations. With this, both Eigen and the compiler can optimize the computations.

The parametric setup also allows for an efficient restriction of a dG or cG function to the Gauss points. Let $A_h\in W_h^{\left(\mathrm{6}\right)}$ and let T∈𝒯_h be again any mesh element with element index $i\in \mathit{\{}\mathrm{1},\mathrm{\dots },N_{\mathrm{el}}\mathit{\}}$. Assume that we want to evaluate A_h in the 3×3 Gauss points ${\widehat{\mathit{\chi }}}_q\in (\mathrm{0},\mathrm{1}{)}^{\mathrm{2}}$; see Sect. 3.5. It holds ${\mathit{\chi }}_q^i:={\mathbf{T}}_T\left({\widehat{\mathit{\chi }}}_q\right)$ and hence

$$\begin{array}{}\text{(25)}& \begin{array}{rl}{A}_h\left({\mathit{\chi }}_q^i\right)& =\sum _{l=\mathrm{1}}^{\mathrm{6}}{\mathbf{A}}_{i,l}{\mathrm{\Psi }}_{\mathrm{l}}\left({\mathbf{T}}_T^{-\mathrm{1}}\right({\mathit{\chi }}_q^i\left)\right)\\ & =\sum _{l=\mathrm{1}}^{\mathrm{6}}{\mathbf{A}}_{i,l}{\widehat{\mathrm{\Psi }}}_{\mathrm{l}}\left({\widehat{\mathit{\chi }}}_q\right)=:{\mathbf{A}}_i^i.\end{array}\end{array}$$

That is, by working with the pulled-back function ${\widehat{\mathrm{\Psi }}}_{\mathrm{l}}$ on the reference element, ${\widehat{\mathrm{\Psi }}}_{\mathrm{l}}$ only needs to be evaluated on the fixed points ${\widehat{\mathit{\chi }}}_q$. Furthermore, by the linearity of the basis representation, the mapping of the local coefficients ${\mathbf{A}}_{i,\mathrm{1}},\mathrm{\dots },{\mathbf{A}}_{i,\mathrm{6}}$ of the dG vector to the values of ${\mathbf{A}}_i^G\in {\mathbb{R}}^{\mathrm{9}}$ in the 9 Gauss points on the element T can be written as a matrix-vector product,

$$\begin{array}{}\text{(26)}& {\mathbf{A}}_i^G={\mathbf{A}}_{i,\cdot }\cdot {\mathbf{G}}_{\mathrm{\Psi }}^{\mathrm{9},\mathrm{6}},[{\mathbf{G}}_{\mathrm{\Psi }}^{\mathrm{9},\mathrm{6}}{]}_{l,q}={\widehat{\mathrm{\Psi }}}_{\mathrm{l}}({\widehat{\mathit{\chi }}}_q),\end{array}$$

with a fixed matrix ${\mathbf{G}}_{\mathrm{\Psi }}^{\mathrm{9},\mathrm{6}}\in {\mathbb{R}}^{\mathrm{9}\times \mathrm{6}}$. The matrices ${\mathbf{G}}_{\mathrm{\Psi }}^{q,s}$ for possible dG degrees with s local unknowns and for supported choices of the Gauss quadrature rule with q points are precomputed and directly inlined into the code to allow for an optimization by Eigen and the compiler. The matrices ${\mathbf{G}}_{\mathrm{\Psi }}^{s,q}$ and similar code are auto-generated by Python scripts to allow for easy extension.

Algorithm 1 Evaluation of the viscosities in the 9 Gauss points in the case of biquadratic velocities.

1:

const Array<9> Ag = A.row(i) * Gpsi<9,3>; // restr. k-th element to GP

2:

const Array<9> Hg = H.row(i) * Gpsi<9,3>; // restr. ice height to GP

3:

const Array<9> E11g = E11.row(i) * Gpsi<9,8>; // restr. strain tensor to GP

4:

const Array<9> E12g = E12.row(i) * Gpsi<9,8>; // restr. strain tensor to GP

5:

const Array<9> E22g = E22.row(i) * Gpsi<9,8>; // restr. strain tensor to GP

6:

7:

const Array<9> zeta = Pstar * Hg * (-C * (1-Ag)).exp() /

8:

(Dmin * Dmin + 1.25 * (E11g+E22g).square() + 1.5 * E11g * E22g + E12g.square()).sqrt();

Another challenge for an efficient implementation is the evaluation of the integrals that are required to determine the viscosity within the VP model (Eq. 4) in the mEVP iteration (see Eq. 19),

$$\begin{array}{rl}\mathit{\zeta }& ={\displaystyle \frac{P^{\star }\cdot H_h\cdot \exp (-C(\mathrm{1}-A_h\left)\right)}{\sqrt{{\mathrm{\Delta }}_{\text{min}}^{\mathrm{2}}+\text{tr}({\mathit{ϵ}}_h{)}^{\mathrm{2}}+\frac{\mathrm{2}}{e^{\mathrm{2}}}{\mathit{ϵ}}_h^{\prime }:{\mathit{ϵ}}_h^{\prime }}}}\\ & ={\displaystyle \frac{P^{\star }\cdot H_h\cdot \exp (-C(\mathrm{1}-A_h\left)\right)}{\sqrt{{\mathrm{\Delta }}_{\text{min}}^{\mathrm{2}}+\frac{\mathrm{5}}{\mathrm{4}}({\mathit{ϵ}}_{\mathrm{11}}+{\mathit{ϵ}}_{\mathrm{22}}{)}^{\mathrm{2}}+\frac{\mathrm{3}}{\mathrm{2}}{\mathit{ϵ}}_{\mathrm{11}}{\mathit{ϵ}}_{\mathrm{22}}+{\mathit{ϵ}}_{\mathrm{12}}^{\mathrm{2}}}}}.\end{array}$$

With $i\in \mathit{\{}\mathrm{1},\mathrm{\dots },N_T\mathit{\}}$ again denoting the element index, the example in Algorithm 1 illustrates the evaluation of the viscosities in the 9 Gauss points in the case of biquadratic velocities, a strain tensor with 8 local unknowns, i.e., $E_h\in [W_h^{\left(\mathrm{8}\right)}{]}^{\mathrm{2}\times \mathrm{2},\mathrm{sym}}$, and tracers discretized as dG(1) functions in $W_h^{\left(\mathrm{3}\right)}$. The above implementation is close to the mathematical notation which simplifies the implementation of model variations. Long expressions such as those in the last line also allow Eigen to vectorize operations efficiently.

4.2 Evaluation of the weak formulations on parametric meshes

A substantial part of the computational effort is due to the mapping of the reference element $\widehat{T}$ onto the mesh elements T∈𝒯_h; compare Sect. 3.1 and 3.5. We discuss the details of an efficient implementation for one specific term in the mEVP momentum Eq. (19), namely the evaluation of $({\mathit{\sigma }}_n^{\left(p\right)},\mathrm{\nabla }{\mathbf{\Phi }}_i{)}_{\mathrm{\Omega }}$, whose discretization has already been given in Eq. (23). In the following, we will omit all indices referring to the time step and the mEVP iteration count.

At the heart of $({\mathit{\sigma }}_n^{\left(p\right)},\mathrm{\nabla }{\mathbf{\Phi }}_i{)}_{\mathrm{\Omega }}$ is the integration of the symmetric stress tensor multiplied with the gradient of the (vector-valued) test function ${\mathbf{\Phi }}_i=({\mathrm{\Phi }}_i^x,{\mathrm{\Phi }}_i^y)$. Pulling this term back from an element T∈𝒯_h onto the reference element $\widehat{T}$, we obtain

$$\begin{array}{}\text{(27)}& \begin{array}{rl}& (\mathit{\sigma },\mathrm{\nabla }{\mathbf{\Phi }}_i{)}_T\\ & =\underset{\widehat{T}}{\int }\text{det}\left(\widehat{\mathrm{\nabla }}{\mathbf{T}}_T\right(\widehat{x}\left)\right)\widehat{\mathrm{\nabla }}{\mathbf{\Phi }}_i\left(\widehat{x}\right)\widehat{\mathrm{\nabla }}{\mathbf{T}}_T^{-T}\left(\widehat{x}\right):\widehat{\mathit{\sigma }}\left(\widehat{x}\right)\text{d}\widehat{x}.\\ & i=\mathrm{1},\mathrm{\dots },N_{\mathrm{loc}}^{\text{cG}}\end{array}\end{array}$$

Here, $N_{\mathrm{loc}}^{\text{cG}}:=(r+\mathrm{1}{)}^{\mathrm{2}}$ is the local number of cG degrees of freedom and $\mathbf{A}:\mathbf{B}:={\sum }_{i,j}{\mathbf{A}}_{ij}{\mathbf{B}}_{ij}$ is the full contraction of rank-2 tensors. Locally on the element T∈𝒯_h, symmetric stresses $\mathit{\sigma }\in [W_h^{\left(s\right)}{]}^{\mathrm{2}\times \mathrm{2},sym}$ and the element map's gradient $\widehat{\mathrm{\nabla }}{\mathbf{T}}_T$ are given in the dG and cG basis as

$$\begin{array}{}\text{(28)}& \begin{array}{rl}\mathit{\sigma }\left(\widehat{x}\right){|}_T& =\sum _{j=\mathrm{1}}^s\underset{={\mathit{\sigma }}_{T,j}}{\underbrace{\left(\begin{array}{cc}{\mathit{\sigma }}_{T,j}^{\mathrm{11}}& {\mathit{\sigma }}_{T,j}^{\mathrm{12}}\\ {\mathit{\sigma }}_{T,j}^{\mathrm{12}}& {\mathit{\sigma }}_{T,j}^{\mathrm{22}}\end{array}\right)}}{\mathrm{\Psi }}_j\left(\widehat{x}\right),\\ \widehat{\mathrm{\nabla }}{\mathbf{T}}_T\left(\widehat{\mathbf{x}}\right)& =\sum _{k=\mathrm{1}}^{\mathrm{4}}\underset{={\mathbf{x}}_{T,k}}{\underbrace{\left(\begin{array}{c}{x}_{T,k}^{\mathrm{1}}\\ x_{T,k}^{\mathrm{2}}\end{array}\right)}}\underset{=\widehat{\mathrm{\nabla }}{\mathrm{\Phi }}_k(\widehat{\mathbf{x}}{)}^T}{\underbrace{\left(\begin{array}{cc}{\partial }_{\widehat{x}}{\mathrm{\Phi }}_k\left(\widehat{\mathbf{x}}\right)& {\partial }_{\widehat{y}}{\mathrm{\Phi }}_k\left(\widehat{\mathbf{x}}\right)\end{array}\right)}}\end{array},\end{array}$$

with the ${\mathbf{x}}_{T,k}=(x_{T,k}^{\mathrm{1}},x_{T,k}^{\mathrm{2}})\in {\mathbb{R}}^{\mathrm{2}}$ being the four corner nodes of the element T. Approximating Eqs. (27) and (28) by Gauss quadrature with n_Q points ${\widehat{\mathbf{x}}}_{\mathrm{1}},\mathrm{\dots },{\widehat{\mathbf{x}}}_{n_Q}\in \widehat{T}$ and weights ${\mathit{\omega }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\omega }}_{n_Q}$ yields

$$\begin{array}{}\text{(29)}& \begin{array}{rl}& (\mathit{\sigma },\mathrm{\nabla }{\mathit{\varphi }}_i{)}_T\\ & \approx \sum _{j=\mathrm{1}}^s\underset{=:{\mathbf{X}}_{i,j}}{\underbrace{\sum _{q=\mathrm{1}}^{n_Q}{\mathit{\omega }}_q\text{det}\left(\widehat{\mathrm{\nabla }}\mathbf{T}\right({\widehat{\mathbf{x}}}_q\left)\right){\mathrm{\Psi }}_j\left({\widehat{\mathbf{x}}}_q\right)\widehat{\mathrm{\nabla }}{\mathrm{\Phi }}_i\left({\widehat{\mathbf{x}}}_q\right)\widehat{\mathrm{\nabla }}\mathbf{T}({\widehat{\mathbf{x}}}_q{)}^{-T}}}:{\mathit{\sigma }}_j.\\ & i=\mathrm{1},\mathrm{\dots },N_{\mathrm{loc}}^{\text{cG}}\end{array}\end{array}$$

The computational effort of the above equation is substantial. The Jacobian $\widehat{\mathrm{\nabla }}{\mathbf{T}}_T$ needs to be assembled $n_Q\cdot s\cdot N_{\mathrm{loc}}^{\text{cG}}$ times (cf. Eq. 29), and its inverse and determinant need to be computed. For the second-order case cG(2) with n_Q=9, $N_{\mathrm{loc}}^{\text{cG}}=\mathrm{9}$, and s=8, more than 15 000 floating point operations are required on each element.

The entries of the 2×2 matrices X_i,j, however, do not depend on the solution but only on the mesh elements T∈𝒯_h. A closer analysis further reveals that ${\mathbf{X}}_{i,j}^{\mathrm{1}}:={\mathbf{X}}_{i,j}^{\mathrm{11}}={\mathbf{X}}_{i,j}^{\mathrm{21}}$ and ${\mathbf{X}}_{i,j}^{\mathrm{2}}:={\mathbf{X}}_{i,j}^{\mathrm{12}}={\mathbf{X}}_{i,j}^{\mathrm{22}}$. Hereby, the complete scalar product with Gauss approximation is evaluated as

$$\begin{array}{}\text{(30)}& \begin{array}{rl}{\left(\mathit{\sigma },\mathrm{\nabla }\left(\begin{array}{c}{\mathrm{\Phi }}_i^x\\ \mathrm{0}\end{array}\right)\right)}_T& \approx {\mathbf{X}}_T^{\mathrm{1}}{\mathit{\sigma }}_T^{\mathrm{11}}+{\mathbf{X}}_T^{\mathrm{2}}{\mathit{\sigma }}_T^{\mathrm{12}},\\ {\left(\mathit{\sigma },\mathrm{\nabla }\left(\begin{array}{c}\mathrm{0}\\ {\mathrm{\Phi }}_i^y\end{array}\right)\right)}_T& \approx {\mathbf{X}}_T^{\mathrm{1}}{\mathit{\sigma }}_T^{\mathrm{12}}+{\mathbf{X}}_T^{\mathrm{2}}{\mathit{\sigma }}_T^{\mathrm{22}},\end{array}\end{array}$$

with matrices ${\mathbf{X}}_T^{\mathrm{1}},{\mathbf{X}}_T^{\mathrm{2}}\in {\mathbb{R}}^{N_{\mathrm{loc}}^{\text{cG}}\times s}$. The computational effort shrinks then to $\mathrm{4}{N}_{\mathrm{loc}}^{\text{cG}}\times s^{\mathrm{2}}$ operations, which in the case of cG(2) amounts to about 2300 operations. The matrices ${\mathbf{X}}_T^{\mathrm{1}}$ and ${\mathbf{X}}_T^{\mathrm{2}}$ can be precomputed and stored for each mesh element. Their small size makes them, furthermore, amenable for efficient caching although additional storage is needed. Section 5.3.3 presents a numerical study on the effective performance of the alternatives, i.e., using precomputed matrices or computation of all terms on the fly.

The same technique can be applied to all further terms of Eq. (17). For some of them the computational savings of precomputing per element terms are even more substantial. This is in particular true if the inverse of the block-diagonal dG mass matrix is required, such as in the mEVP iteration (19).

Figure 3 Parallel processing of the vectors on a small mesh with 12 elements, N_x=4 and N_y=3. All blocks of one color can be processed in parallel without memory conflicts. From left to right: node-wise operations, local element-wise operations, edge-wise operations on e^x edges (blue) and e^y edges (green), and operations writing on biquadratic cG(2) vectors.

4.3 OpenMP parallelization

In each MPI task, only topologically regular rectangular meshes are considered that consist of $N_{\mathrm{el}}:=N_x\times N_y$ elements. As the complete numerical workflow is based on explicit integrators, OpenMP parallelization is easily realized. Depending on the specific task, a different coloring of the mesh elements (or mesh edges) is utilized to avoid any memory conflicts.

Node-wise.

Vector operations (such as sums and entry-wise products) are parallel with respect to the major index referring to the node.

Element-wise.

Operations such as the stress update in Eq. (19) within the mEVP iteration in Eq. (17) are parallel with respect to the mesh element. This also includes the element-wise terms (H_hv,∇ψ)_T of the transport scheme in Eq. (16) where no communication is involved and also the projection of the strain rate tensor from the cG to the dG space:

$$(J_T{\mathbf{E}}_h,\mathrm{\Psi }{)}_{\widehat{T}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(J_T\left(\widehat{\mathrm{\nabla }}\widehat{\mathbf{v}}\right[\widehat{\mathrm{\nabla }}{\mathbf{T}}_T{]}^{-\mathrm{1}}+\left[\widehat{\mathrm{\nabla }}{\mathbf{T}}_T{]}^{-T}\widehat{\mathrm{\nabla }}{\widehat{\mathbf{v}}}^T\right),\mathrm{\Psi }{)}_{\widehat{T}}.$$

Edge-wise.

The edge integrals in Eq. (16) are assembled in two sweeps. First, all horizontal edges e_x∈𝒯_h are computed,

$$\sum _{i_x=\mathrm{1}}^{N_x}\sum _{i_y=\mathrm{0}}^{N_y}{〈\mathit{\left\{}\mathit{\right\{}{H}_h\mathit{\left\}}\mathit{\right\}},\mathbf{v}\cdot {\mathit{n}}_e\left[\right[\mathit{\psi }\left]\right]〉}_{e_{i_x,i_y}^x}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{〈|\mathbf{v}\cdot {\mathit{n}}_e|\cdot \left[\right[H_h\left]\right],\left[\right[\mathit{\psi }\left]\right]〉}_{e_{i_x,i_y}^x},$$

and the outer (in the x direction) is run in parallel as the integral on an edge $e_{i_x,i_y}^x$ will affect the two elements atop and below it. Then, a second sweep, parallelized in the y direction, performs the computation for the e^y edges.

When updating cG vectors, e.g., in the stress update (cf. Eq. 23), more care is required. We use a row-wise coloring of the elements and perform the update in two sweeps. Figure 3 summarizes the parallel processing of the mesh.

Remark 2 (towards GPU acceleration).

Our finite-element discretization requires a large number of per element computations with a substantial amount of computations for each one. Furthermore, the computational costs increase substantially with the order; see Sect. 4.2. Only local coupling between adjacent elements thereby exists since an explicit time stepping and mEVP iterations are used. This makes the problem well suited for a GPU parallelization where thousands of independent computations are required to fully utilize a state-of-the-art GPU and even more when multiple GPUs are combined in a node. The current implementation has already been designed with a GPU implementation in mind. Its realization is planned as a next step.

5 Numerical experiments

In this section we will present a set of experiments to validate our discretization. We will thereby first only study the accuracy of the advection before considering the full mEVP scheme.

Figure 4 Advection test case I: convergence rates on uniform (red) and distorted (blue) meshes.

5.1 Validating the higher-order transport scheme

5.1.1 Advection test case I: transport of an initially smooth bump

On the domain $\mathrm{\Omega }=(\mathrm{0},L_x)\times (\mathrm{0},L_y)$ with L_x=409 600 and L_y=512 000 we advect the initially smooth bump,

$$\begin{array}{rl}& H_{\text{in}}\left(\mathbf{x}\right)=\left\{\begin{array}{ll}\exp (-\frac{\mathrm{1}}{\mathrm{1}-r\left(\mathbf{x}\right)})& r\left(\mathbf{x}\right)<\mathrm{1}\\ \mathrm{0}& r\left(\mathbf{x}\right)\ge \mathrm{0}\end{array}\right.,\\ & r\left(\mathbf{x}\right)=\mathrm{40}\Vert {\displaystyle \frac{\mathbf{x}}{L_x}}-({\displaystyle \frac{\mathrm{1}}{\mathrm{4}}},{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{)}^T{\Vert }^{\mathrm{2}},\end{array}$$

with the stationary, rotational velocity field,

$$\mathbf{v}\left(\mathbf{x}\right)={\displaystyle \frac{\mathit{\pi }}{L_x}}\left(\begin{array}{c}\mathrm{2}{\mathbf{x}}_{\mathrm{2}}-L_x\\ L_x-\mathrm{2}{\mathbf{x}}_{\mathrm{1}}\end{array}\right).$$

The problem is run in the time interval $T=[\mathrm{0},L_x]$ such that one complete revolution of the bump is performed. We compute the test case on a sequence of meshes consisting of $N_x^{\left(l\right)}\times N_y^{\left(l\right)}$ elements and $N_T^{\left(l\right)}$ time steps using

$$\begin{array}{rl}& N_x^{\left(l\right)}=\mathrm{24}\cdot {\mathrm{2}}^{l-\mathrm{1}},\\ & N_y^{\left(l\right)}=\mathrm{26}\cdot {\mathrm{2}}^{l-\mathrm{1}},\\ & N_T^{\left(l\right)}=\mathrm{200}\cdot {\mathrm{2}}^{l-\mathrm{1}}\cdot (r+\mathrm{1}{)}^{-\mathrm{2}},l=\mathrm{1},\mathrm{2},\mathrm{\dots },\end{array}$$

where $r\in \mathit{\{}\mathrm{0},\mathrm{1},\mathrm{2}\mathit{\}}$ is the degree of the dG(r) approach. The coarsest discretization consists of $\mathrm{24}\cdot \mathrm{26}=\mathrm{624}$ elements of approximate size 17 km×19.6 km each and a time step Δt=512 s. This results in a CFL constant lower than $\mathrm{0.5}\cdot (r+\mathrm{1}{)}^{-\mathrm{2}}$, which is sufficient for a robust discretization. Next to these uniform rectangular meshes, we use a sequence of distorted meshes to model the effect one encounters in a mesh parametrization of the sphere; see Fig. 5. The nodes x_i,j are in this case given by

$$\begin{array}{rl}& {\mathbf{x}}_{i,j}=\left(\begin{array}{c}\frac{i\cdot L_x}{N_x}+\frac{\mathrm{1}}{\mathrm{20}}\sin \left(\frac{i\cdot \mathrm{3}\mathit{\pi }}{N_x}\right)\sin \left(\frac{j\cdot \mathit{\pi }}{N_y}\right)\\ \frac{j\cdot L_y}{N_y}+\frac{\mathrm{1}}{\mathrm{20}}\sin \left(\frac{i\cdot \mathrm{2}\mathit{\pi }}{N_x}\right)\sin \left(\frac{j\cdot \mathrm{2}\mathit{\pi }}{N_y}\right)\end{array}\right),\\ & \text{ for }i=\mathrm{0},\mathrm{\dots },N_x,\text{ and }j=\mathrm{0},\mathrm{\dots },N_y.\end{array}$$

Through the periodicity of the domain, the exact solution at time T=L_x equals the initial condition. We measure the scaled L² error by

$$\text{err}={\displaystyle \frac{\mathrm{1}}{L_x}}\Vert H_{h,\mathrm{\Delta }t}\left(T\right)-H_{\text{in}}{\Vert }_{L^{\mathrm{2}}\left(\mathrm{\Omega }\right)}.$$

The scaling factor $\mathrm{1}/L_x$ accounts for the drift-error accumulation that is expected to be dependent on the length of the advection in space.

Figure 4 shows the convergence behavior for the different meshes and degrees r. We observe the expected convergence rate of $\mathcal{O}\left(\right|\mathrm{\Delta }x{|}^{\frac{\mathrm{1}}{\mathrm{2}}})$ for dG(0) (cf.  Di Pietro and Ern, 2012, Theorem 3.7) and even super-convergent second order instead of $\mathcal{O}\left(\right|\mathrm{\Delta }x{|}^{\mathrm{1}+\frac{\mathrm{1}}{\mathrm{2}}})$ for dG(1) and super-convergent third order instead of $\mathcal{O}\left(\right|\mathrm{\Delta }x{|}^{\mathrm{2}+\frac{\mathrm{1}}{\mathrm{2}}})$ for dG(2) (see  Di Pietro and Ern, 2012, Theorem 3.13). Distortion of the meshes slightly increases the error constant, but the convergence order is not affected, as expected.

Figure 5 Advection test case I: visualization of the coarse meshes and the initial dG(0) solution. (a) Regular rectangular mesh. (b) Distorted parametric mesh.

5.1.2 Advection test case II: transport in a circular annulus

The domain of the second test case is a circular annulus with inner radius r₀=100 km and outer radius r₁=250 km; see Fig. 6. The parametric mesh is constructed by mapping a uniform rectangular mesh onto the ring using the map

$$\mathbf{T}\left({\mathbf{x}}_i\right):=(r_{\mathrm{0}}+(r_{\mathrm{1}}-r_{\mathrm{0}}\left){\displaystyle \frac{i_y}{N_y}}\right)\left(\begin{array}{c}\cos \left(\frac{\mathrm{2}\mathit{\pi }\cdot i_x}{N_x}\right)\\ -\sin \left(\frac{\mathrm{2}\mathit{\pi }\cdot i_x}{N_x}\right)\end{array}\right).$$

The divergence free stationary velocity field for the transport is given by

$$\mathbf{v}\left(\mathbf{x}\right)={\displaystyle \frac{\mathrm{2}\mathit{\pi }\mathrm{m}}{\mathrm{250}\mathrm{000}\mathrm{s}}}\cdot \left(\begin{array}{c}{\mathbf{x}}_{\mathrm{2}}\\ -{\mathbf{x}}_{\mathrm{1}}\end{array}\right),$$

and it moves the initial conditions uniformly along the domain.

Figure 6 Advection test case II: computational domain and mesh (a) for the transport in a ring. The boundary lines marked in red are the periodic boundaries of the regular reference mesh. On panel (b) we show the initial solution.

One complete revolution around the annulus is achieved in T=250 000 s. The initial field consists of four objects with different regularity: a smooth C^∞ bump centered at $(-\mathrm{175}\mathrm{km},\mathrm{0}\mathrm{km})$ of radius 50 km (on the left), a continuous C⁰ pyramid centered at (175 km,0 km) with radius 50 km (on the right), and two discontinuous discs with radius 50 km at $(\mathrm{0}\mathrm{km},-\mathrm{175}\mathrm{km})$ (on the bottom) and (0 km,175 km) (on the top). The last one has a Pac-Man-shaped omission.

Figure 7 Difference (absolute value) between the initial condition and the solution after one revolution for dG(0) (a), dG(1) (b), and dG(2) (c).

Figure 7 shows the absolute error between the initial condition and the solution at time T. The results verify the superiority of higher-order methods with respect to numerical diffusion. While the finite-volume dG(0) discretization is highly diffusive and inaccurate, both dG(1) and dG(2) preserve all shapes well with very small errors for the smooth ones and small but sharp interface errors for the disc and the Pac-Man-like shape.

5.2 Validating the mEVP framework

To validate the complete dynamical mEVP framework and to assess the higher-order discretization, we study the VP benchmark problem that has been introduced in Mehlmann and Richter (2017) and investigated in Mehlmann et al. (2021b) to compare different sea ice realizations with respect to their ability to depict linear kinematic features.

Figure 8 Setup of the benchmark problem. Forcing by means of a steady circular ocean current (in blue) and a divergent wind field (in red) that is moving diagonally across the domain (with side length 512 km). Initially the ice has an average height of 0.3 m (with a periodic variation) and concentration A=1. The simulation is run for a period of 2 d, where the wind's center travels from the midpoint to the upper right (red arrow).

Table 1 Default values of the VP model used to define the benchmark (Mehlmann and Richter, 2017) as well as default mEVP parameters used in all numerical test cases.

The setup of the benchmark is described in Fig. 8. On the domain of size 512 km×512 km, the ice has an initial concentration of one and an average height of 0.3 m with small amplitude oscillations of 0.005 m and with a wavelength of 105 km in the x direction and 210 km in the y direction; i.e., we have

$$\begin{array}{rl}\mathbf{v}{|}_{t=\mathrm{0}}& =\mathrm{0},A\left(\mathbf{x}\right){|}_{t=\mathrm{0}}=\mathrm{1},\\ H\left(\mathbf{x}\right){|}_{t=\mathrm{0}}& =\mathrm{0.3}\mathrm{m}+\mathrm{0.005}\mathrm{m}\left(\sin \right({\displaystyle \frac{\mathrm{6}{\mathbf{x}}_{\mathrm{1}}}{\mathrm{100}\mathrm{km}}})+\sin ({\displaystyle \frac{\mathrm{3}{\mathbf{x}}_{\mathrm{2}}}{\mathrm{100}\mathrm{km}}}\left)\right).\end{array}$$

The forcing in the benchmark problem consists of a rotational ocean forcing,

$${\mathbf{v}}_{\mathrm{o}}\left(\mathbf{x}\right)={\displaystyle \frac{\mathrm{0.01}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}}{\mathrm{512}\mathrm{km}}}\left(\begin{array}{c}\mathrm{2}{\mathbf{x}}_{\mathrm{2}}-\mathrm{512}\mathrm{km}\\ \mathrm{512}\mathrm{km}-\mathrm{2}{\mathbf{x}}_{\mathrm{1}}\end{array}\right),$$

and a rotational divergent wind field with a center m(t) that is moving along the diagonal of the domain:

$$\begin{array}{rl}& {\mathbf{v}}_{\mathrm{a}}\left(\mathbf{x}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{100}}}\exp \left(\mathrm{1}-{\displaystyle \frac{|\mathbf{x}-\mathbf{m}(t){|}_{\mathrm{2}}}{\mathrm{100}\mathrm{km}}}\right)\left(\begin{array}{cc}\cos \left(\mathit{\alpha }\right)& \sin \left(\mathit{\alpha }\right)\\ -\sin \left(\mathit{\alpha }\right)& \cos \left(\mathit{\alpha }\right)\end{array}\right)\\ & \cdot (\mathbf{x}-\mathbf{m}(t\left)\right),\mathbf{m}\left(t\right)=(\mathrm{256}\mathrm{km}+t\cdot \mathrm{51.2}\mathrm{km}{\mathrm{d}}^{-\mathrm{1}})\left(\begin{array}{l}\mathrm{1}\\ \mathrm{1}\end{array}\right).\end{array}$$

The different parameters of the VP model and the mEVP iteration are given in Table 1.

Figure 9 Shear deformation (log ₁₀) for the benchmark at time T=2 d. We compare different combinations of cG velocity discretizations (linear and quadratic) with dG advection discretization (constant, linear, and quadratic). All computations are run on meshes with the grid spacing h=2 km.

Figure 9 shows the shear $S\left(\mathit{ϵ}\right):=\sqrt{({\mathit{ϵ}}_{\mathrm{11}}-{\mathit{ϵ}}_{\mathrm{22}}{)}^{\mathrm{2}}+\mathrm{4}{\mathit{ϵ}}_{\mathrm{12}}^{\mathrm{2}}}$ at final time T=2 d on a mesh with h=2 km spacing. Results are presented for all combinations of the velocity discretization (low-order cG(1) and high-order cG(2)) and advection schemes dG(0), dG(1), and dG(2). The results are in very good agreement with earlier results for the benchmark (Mehlmann and Richter, 2017) and also with recent numerical studies that consider some of the most widely used sea ice models (Mehlmann et al., 2021b).

For all combinations of velocity and dG degrees, linear kinematic features are well resolved and the deformation field is stable (for a detailed discussion, we refer to Sect. 5.2.1). The results of Fig. 9 suggest that the role of the advection scheme is minor in comparison to the discretization of the velocities. This will be discussed in Sect. 5.2.2. While the cG(1) results in the top row of Fig. 9 are comparable to the B-grid staggerings (cG(1) is the finite-element equivalent of this) given in Mehlmann et al. (2021b, Fig. 6), the high-order cG(2) results show patterns that are at least as resolved as the CD-grid results in Mehlmann et al. (2021b, Figs. 6 and 7).

Figure 10 Total length and number of linear kinematic features detected for the low-order cG(1)-dG(0) and high-order cG(2)-dG(2) schemes. For comparison, we show selected simulation results on quadrilateral meshes from various models taken from Mehlmann et al. (2021b).

5.2.1 Resolution of LKFs

To further investigate the ability of the cG-dG framework to resolve linear kinematic features, we use code provided by Hutter et al. (2019)¹ that identifies LKFs from the shear deformation rate field. The original scripts have been slightly modified in the following manner: the resolution of the uniform quadrilateral mesh onto which the outputs are initially projected has been increased from 256×256 to 512×512 to account for the fact that the higher-order dG discretizations carry subgrid information that would otherwise be lost. The length of LKFs, which the scripts measure in pixels, was adjusted accordingly. Figure 10 shows the results for a selection of the originally published data sets (Mehlmann et al., 2021b, a) together with the results of the low-order cG(1)-dG(0) and high-order cG(2)-dG(2) simulations performed with the proposed discretization. The low-order results are consistent with the data published by Mehlmann et al. (2021b). In particular, the results agree with those obtained with Gascoigne (B-grid) (Braack et al., 2021), which is based on the same discretization. The high-order cG(2) scheme of our discretization can resolve substantially more (and longer) features on coarser meshes, demonstrating the advantage of higher-order schemes.

We note that the exact number and length of features is affected by the chosen mEVP parameters. We have not applied any fine-tuning here but use the values given in Table 1 for all meshes and all discretizations. Moreover, a direct comparison of results obtained on quadrilaterals and triangles is difficult, so we refrain from a more detailed analysis. We refer to Mehlmann et al. (2021b, Sect. 6) for an in-depth discussion of these aspects.

Table 2 Computational time for the benchmark test case using the 4 km mesh. The specified times include the advection and the momentum equation for the complete simulation spanning from t=0 to t=2 d. Times for input/output are not included. Thirty-two parallel threads are used.

Table 2 indicates computational times for simulating the benchmark problem on the 4 km mesh for all different cG/dG combinations for velocity and tracers. The simulations are run on an AMD EPYC 7662 64-core processor at 3.20 GHz using 32 cores. The given times are the total times for the complete benchmark run, and they include everything apart from the input/output operations. The impact of the degree of the tracers is negligible as the advection scheme only takes a little fraction of the overall runtime. Scaling from linear to quadratic velocities is less severe than expected. This is due to the larger fraction of local work at high order and associated with this a better possibility of vectorization by Eigen as well as better efficiency of parallelization.

5.2.2 The role of the advection scheme

In the literature, the role of the advection scheme in the VP model and in particular for the LKF resolution is not entirely clear (again, see Mehlmann et al., 2021b, Sect. 6). We therefore compare for the above benchmark results for a high-order momentum discretization (biquadratic velocities) and dG(0), dG(1), and dG(2) advection to shed further light on the effect of the advection.

Figure 11 Visualization of the ice concentration on meshes with a spacing of h=8 km for dG(0) (finite volumes), dG(1), and dG(2).

Figure 11 shows the ice concentration at time T=2 d for the benchmark problem run on a mesh with spacing h=8 km. The velocity is discretized bi-quadratically with cG(2), and the tracers are represented as discontinuous dG(0), dG(1), and dG(2) functions. The elevated surfaces in the middle and on the right in Fig. 11 show the additional information that is gained by higher-order approaches. The finite-volume dG(0) discretization only gives average values in each element, while the dG(1) discretization includes the slope of the tracers, and starting with dG(2) further information on the subgrid-scale, e.g., on the curvature, is also represented. The effect of the dG degree on the representation of the sea ice drift (and derived values like the shear deformation) is less drastic, as shown in Fig. 9. In comparison to the choice of the velocity discretization, the effect of the tracer discretization on the representation of LKFs is small. This is also the conclusion of Mehlmann et al. (2021b). Figure 12 shows the results of the LKF detection code by Hutter et al. (2019) for different advection schemes. While the low-order discretization with cG(1) velocities is hardly affected by the advection discretization, the impact on the cG(2) high-order scheme is larger. Here, the lowest-order upwind scheme dG(0) yields the most and longest LKFs.

Figure 12 Effect of the advection scheme on resolving linear kinematic features. (a) Total length of detected LKFs. (b) Number of LKFs.

5.2.3 Stability of the mixed finite-element formulation

Section 3.4 introduced the mEVP iteration as a mixed finite-element formulation, and in particular Remark 3.4 discussed the optimal choice of the velocity and stress spaces. In this section, we demonstrate the effect of this choice on the results. We again consider the benchmark problem of Sect. 5.2 using a mesh with h=2 km. The tracers are discretized with dG(1); all parameter values are as given in Table 1.

Figure 13 Shear at final time T=2 for different choices of the velocity space $V_h^{\left(r\right)}$ and stress space $[W_h^{\left(s\right)}{]}^{\mathrm{2}\times \mathrm{2},\mathrm{sym}}$ (detail from the simulation domain). (a, b) Unstable finite-element pair with stresses one degree below the velocity space. (c, d) Stable combinations.

Figure 13 shows a snapshot of the shear deformation rate at time T=2 d for all different combinations of velocity and stress spaces. The plots clearly show the need to use large-enough stress spaces. For the low-order velocity $V_h^{\left(\mathrm{1}\right)}$ the results based on piece-wise constant stresses $\mathit{\sigma }\in [W_h^{\left(\mathrm{1}\right)}{]}^{\mathrm{2}\times \mathrm{2},\mathrm{sym}}$ (upper left) appears to give reasonable results, in particular when compared to the highly diffusive combination $V_h^{\left(\mathrm{1}\right)}-W_h^{\left(\mathrm{3}\right)}$ (lower left). $V_h^{\left(\mathrm{1}\right)}-W_h^{\left(\mathrm{1}\right)}$, however, is unstable and does not satisfy Eq. (22). It shows oscillations on the level of the mesh elements, while the combination $V_h^{\left(\mathrm{1}\right)}-W_h^{\left(\mathrm{3}\right)}$ is stable. Stress spaces that are too small do not provide sufficient control of the term (σ(v),Ψ) in Eq. (19) and lead to oscillatory stresses; see the upper right plot in Fig. 13. We refer to the discussion in Sect. 3.4.

5.3 Computational efficiency

5.3.1 Shared-memory node-level parallelization

The complete code is parallelized using OpenMP to benefit from shared-memory parallelism available on individual nodes. To evaluate the parallel performance of the code, we run a strong scalability test. The benchmark problem described in Sect. 5.2.1 is run for T=1 h on a mesh with h=2 km spacing which amounts to $N={\mathrm{256}}^{\mathrm{2}}=\mathrm{65}\mathrm{536}$ elements. Using a time step size of Δt=60 s a total of 60 advection time steps and mEVP subiterations with 100 steps each are computed. We discretize the velocity with quadratic cG(2) elements and use dG(2) with 6 local unknowns for ice height and sea ice concentration.

The simulation is run on an AMD EPYC 7662 64-core processor at 3.20 GHz. Figure 14 shows the strong scalability results. The overall runtime drops from 245 s on 1 core to 6.5 s on 64 cores. The parallel efficiency stays very high at about 0.9 when run on up to 16 cores and then slightly drops. The parallel efficiency of the advection scheme is not as good as the efficiency of the mEVP iteration. In this benchmark problem, this is not significant since only two tracers are advected and since there are 100 substeps of the mEVP solver in each advection step. For more complex thermodynamics, the situation will be different and further optimizations appear to be necessary. However, the parallelization effort so far has been restricted to enabling OpenMP. A GPU implementation also provides significant opportunities for better parallel scaling.

Figure 14 Scaling (a) and parallel efficiency (b) of the OpenMP parallelization for the benchmark test case on the 2 km mesh. The mesh has 65 536 elements, and a total of 60 advection steps with 100 momentum substeps within the mEVP iteration are covered by the total time.

5.3.2 Vectorization

Our implementation described in Sect. 4 benefits from the vectorization capabilities of Eigen (Guennebaud et al., 2010). Eigen can in particular exploit the additional computations that arise in higher-order discontinuous Galerkin methods as the dG space $W_h^{\left(s\right)}$ naturally leads to many s×s matrices. Furthermore, the use of parametric meshes makes numerical quadrature with Gauss rules necessary and leads to operations involving n_q×n_cG matrices, where n_q is the number of Gauss points in an element (usually 9) and n_cG is the local number of cG basis functions (4 for cG(1) and 9 for cG(2)). Table 3 shows the effect of CPU-level vectorization.²

Table 3 The effect of vectorization on the computational times for different discretizations of the momentum and the advection equations.

Table 4 Comparison of different mesh structures. Evaluation of the benchmark problem with cG(2)-dG(2) discretization on meshes of size h=8 km, h=4 km, and h=2 km. The benchmark is simulated for T=1 h. Parts (a) and (b) refer to the parametric meshes, where the variational formulation must be evaluated using Gauss quadrature; part (c) corresponds to the uniform mesh, where the exact evaluation of the variational form is hard-coded. Part (b) in comparison to part (a) precomputes and stores the matrices required for the parametric mapping.

The effect is significant and purely based on the design principle in our implementation to use Eigen as much as possible. Computations are run using 16 cores of the AMD EPYC 7662 64-core CPU running at 3.20 GHz, and reported is the average of three consecutive calls.

5.3.3 Computational overhead of parametric meshes

As detailed in Sect. 4, with a parametric mesh the variational formulation needs to be mapped back onto the reference element for numerical quadrature. This has substantial computational overhead compared to fully uniform grids where all essential quantities can be precomputed. However, we described in Sect. 4 that also in the parametric case substantial precomputations are possible.

Table 4 shows memory usage and computation times as a function of the mesh type and for three successively refined meshes. We also provide computational time and memory consumption per element of the mesh. We again solve the benchmark problem of Sect. 5.2.1 but only in the short interval [0,1 h] using 16 CPU cores and a high-order discretization, i.e., cG(2) with dG(2) advection. The results in Table 4 clearly show the superiority of the uniform mesh approach both in terms of computational time and memory consumption. It shows, however, that parametric meshes can be realized either with a comparable memory footprint but with substantial computational overhead or with comparable computational efficiency but with increased memory requirements.

The trade-off that is chosen in an implementation will likely depend on the hardware that is targeted. On common multi-core CPUs with a moderate number of cores, the parametric approach with precomputed matrices seems to be superior. When using many-core systems or moving to GPUs or TPUs, the balance between compute-memory performance may differ. This remains a task for further research in the future.

6 Conclusion

We presented the numerics and implementation of the neXtSIM-DG dynamical core, a new discretization of sea ice dynamics aimed at Earth system models. A key feature is the use of a higher order in terms of local discretization and the consistent use of efficient data structures and modern programming paradigms.

The new framework has been validated in the context of Hibler's established viscous–plastic sea ice model, but the dynamical core is flexible and can accommodate different rheology models.

All advection equations are discretized using discontinuous Galerkin methods. Currently, the methods dG(0), dG(1), and dG(2) have been implemented and validated, but the flexible software concept based on code-generation and pre-assembled matrices for efficient implementation of the variational formulations easily allows an extension to even higher orders. This could become relevant in connection with alternative sea ice rheologies; see, e.g., Dansereau et al. (2016) and Rampal et al. (2016). The momentum equation is discretized using second-order continuous finite elements.

We validated the advection discretization and showed that the theoretically expected orders of convergence are realized in practice. Thereby, by using parametric grids, we achieve great flexibility on the spatial discretization. On the other hand, the underlying structured grid topology allows for an efficient numerical implementation. The momentum equation with an mEVP approximation of the visco-plastic model and was tested using an established benchmark problem (Mehlmann and Richter, 2017). In particular, we showed that the high-order discretization can resolve more LKFs than the established models (Mehlmann et al., 2021b), albeit with a larger number of degrees of freedom.

We also described the implementation of the dynamical core on a shared-memory compute node and parallelization using OpenMP. In our future work, we will add coarser parallelization on distributed clusters with MPI and also parallelization on GPUs.

While MPI parallelization is standard and can be easily accomplished by using topologically simple, structured rectangular grids, we enter new territory with GPU parallelization. However, this is already prepared by using a structured memory design. Moreover, the current implementation of neXtSIM-DG allows an easy choice between a pre-assembly of matrices and an on-the-fly computation of all sizes, which could be beneficial on GPUs to reduce memory bandwidth; see Sect. 5.3.3.