Finite-Element Sea Ice Model (FESIM), version 2

The Finite-Element Sea Ice Model (FESIM), used as a component of the Finite-Element Sea ice Ocean Model, is presented. Version 2 includes the elastic-viscous-plastic (EVP) and viscous-plastic (VP) solvers and employs a flux corrected transport algorithm to advect the ice and snow mean thicknesses and concentration. The EVP part also includes a modified approach proposed recently by Bouillon et al. (2013), which is characterized by an improved stability compared to the standard EVP approach. The model is formulated on unstructured triangular meshes. It assumes a collocated placement of ice velocities, mean thicknesses and concentration at mesh vertices, and relies on piecewise-linear (P1) continuous elements. Simple tests for the modified EVP and VP solvers are presented to show that they may produce very close results provided the number of iterations is sufficiently high.


Introduction
The Finite-Element Sea Ice Model (FESIM) was developed as a component of the Finite-Element Sea Ice Ocean circulation Model (FESOM) (for a recent description see Wang et al., 2014) in 2003.Its basis was the standard zero-layer thermodynamical component, and an elastic-viscous-plastic (EVP) solver coded following Hunke and Dukowicz (1997) and the early version of CICE documentation (see Hunke and Lipscomb, 2008 for the current one).It was the first unstructured-mesh sea ice model used for global oceansea ice simulations.The description of the first version was only available as an internal technical report (Danilov and Iakovlev, 2003, unpublished manuscript) and in a brief form was presented by Timmermann et al. (2009).The P 1 − P 1 (linear polynomials on triangles for velocities and scalars) continuous representation used in the dynamical core led to a very compact code relying on the numerical infrastructure of FESOM.The components of stresses and strain rate tensors are elementwise constant, which makes the numerical implementation very straightforward.
Version 2 of the model is augmented by a new viscousplastic (VP) solver, while the Galerkin least-squares stabilized advection scheme inherited from early versions of FE-SOM is replaced by the FE (finite-element) flux corrected transport (FCT) scheme by Löhner et al. (1987), which warrants better numerical stability.It also contains the new EVP solver by Bouillon et al. (2013), which puts the EVP and VP approaches on the same footing.The model reached a high level of maturity and shows a robust behavior in numerous simulations performed with FESOM (see, e.g., Sidorenko et al., 2011;Wang et al., 2012;Wekerle et al., 2013;Timmermann and Hellmer, 2013;Wang et al., 2014;Sidorenko et al., 2015).It may serve as a prototype for other groups developing unstructured-mesh models intended for large-scale ocean sea ice simulations.
The intention of this paper is to present the description of the dynamical part of the model (momentum balance and tracer advection), and illustrate the performance of the solver algorithms implemented in the model.The thermodynamical part will not be described here, as its implementation is standard (pointwise) and is not affected by the unstructured character of the surface mesh.It follows Parkinson and Washington (1979) and includes a prognostic snow layer (Owens and Lemke, 1990).
Several approaches to sea ice modeling on unstructured meshes have been proposed recently.Hutchings et al. (2004) describe an approach based on a finite-volume (FV) cell-centered discretization.Another finite-volume implementation is that by FVCOM (the unstructured grid Finite Volume Community Ocean Model), which follows CICE (see Hunke and Lipscomb, 2008) but employs cell-vertex discretization, i.e., velocities are on cells (triangles), and tracers are on vertices (see Gao et al., 2011).Next to FESIM, another FE model has been proposed by Lietaer et al. (2008).It relies on linear non-conforming elements for velocities (full velocity vectors are associated with the edges of the triangular mesh) and elementwise constant tracers.We comment on these discretizations later.
Sections 2 and 3 introduce the basic equations and present the description of the model's numerical part.We discuss some aspects of model performance in Sect. 4 and conclude the presentation in Sect. 5.

Governing equations
The sea ice motion equation is Here m is the ice plus snow mass per unit area, C d the ice-ocean drag coefficient, ρ o the water density, a the sea ice concentration, u = (u, v) and u o the ice and ocean velocities, τ the wind stress, H the sea surface elevation, g the acceleration due to gravity and F j = ∂ i σ ij is the force from stresses within the ice.We use Cartesian coordinates for brevity, with i, j = 1, 2 implying x and y directions; the implementation of spherical coordinates will be discussed later.Summation over repeating coordinate indices is implied.The total mass m is with ρ ice and ρ s , respectively, the densities of ice and snow and h ice and h s their mean thicknesses (volumes per unit area).
The internal ice stresses are computed assuming the VP rheology (Hibler, 1979).One writes where is the strain rate tensor, η and ζ are the moduli ("viscosities") and P is the ice strength.Both the stress and the strain rate tensors are symmetric, so they are characterized by only three independent components.The standard VP rheology adopts the following scheme of computing the ice strength P and moduli η and ζ : where e = 2 (the ellipticity parameter) and C = 20; the default values in FESOM for min and p * are min = 2 × 10 −9 s −1 and p * = 27 500 N m −2 .In this scheme, min serves for a viscous regularization of plastic behavior in areas where is very small.The ice strength can be modified as P = P 0 /( + min ) for stresses to remain on the elliptic yield curve even if is small, and we will follow this variant below.We note that multi-category ice implementations (such as CICE; see Hunke and Lipscomb, 2008) use different parameterizations for P 0 , which take into account the distribution of ice over thickness categories.This does not change the basic Eqs. ( 1) and (3).
In our case we deal with three tracers, the concentration a, ice mean thickness (volume per unit area) h ice and snow mean thickness h s .They are advected by the ice velocities and modified through thermodynamical forcing: with S a and S ice the sources related to sea ice melting and freezing, and S s the sources due to snow precipitation and melting.The system (1), ( 3) and ( 7), augmented with an appropriate model of sources and boundary conditions, defines the sea ice model.We use the no-slip boundary conditions for momentum and no-flux condition for tracers at lateral walls.

VP and EVP methods
The well known difficulty in solving the ice momentum equation is related to the internal stress term, which makes this equation very stiff and would require time steps of fractions of a second if stepped explicitly.There are two common ways of handling this difficulty.The first one treats a part of stress divergence in an implicit way, with linearization for the moduli, as suggested by Zhang and Hibler (1997).As mentioned by Lemieux and Tremblay (2009), it does not warrant full convergence, and a full nonlinear solver (for example, a Jacobian-free Newton-Krylov solver; see Lemieux et al., 2012) has to be used for that.This strategy is still too expensive computationally, so the VP solver adopted by us is similar in spirit to that of Zhang and Hibler (1997) (see Sect. 3.4).
The second way is to reformulate the VP approach by adding pseudo-elasticity, which leads to the so-called EVP method.
It raises the order of the system (1) and (3) with respect to time, which makes the CFL (Courant-Friedrichs-Lewy) limitation on the explicit time step much less severe than in the original VP framework.The EVP approach, as proposed by Hunke and Dukowicz (1997) (see also Hunke and Lipscomb, 2008), is described as follows.One first defines the combinations and similar combinations for the strain rate components: In this notation, the EVP approach is where T is the relaxation time.It determines the timescale of transition from elastic behavior to the VP rheology.The default value is T = t/3, where t is the external time step (set by the ocean model).It can be easily seen that the EVP "rheology" becomes equivalent to the VP rheology if the contribution from the time derivatives are negligible on the timescale given by t.The equations for stresses are time stepped together with the momentum Eq. ( 1) at a shorter time step t EVP , so that N EVP = t/ t EVP is a large number (about 100 or more).A caveat of this approach is that by the end of the external time step the stresses may still differ from the VP solution, and the difference may accumulate with time.So, in practice the EVP solution may slightly deviate from the VP one.Because of purely explicit time stepping for the stress-velocity pair (velocity is considered known in stress computations and vice versa), the EVP approach must respect the CFL limitation on the subcycling time step t EVP (see Hunke and Dukowicz, 1997;Hunke, 2001).It can be circumvented by limiting "viscosities" (ζ = P 0 /2( + min ), η = ζ /e 2 ) so that they stay below some level (see Hunke, 2001) , where C lim is the limiting constant and x the grid cell size.However, on unstructured meshes this can modify solutions simply because of varying resolution (see the discussion by Losch and Danilov, 2012).Limiting is therefore not used by us.The stability condition then demands that t EVP remains small.Note that the limitation on t EVP becomes more restrictive for finer meshes, and would require to use a larger N EVP .
If not observed, the CFL limitation may lead to noisy fields of velocity divergence and viscosities in practical applications in the areas where is low.The code remains stable in most cases (because of stability added through time stepping, see further) and produces relatively smooth results for the ice thickness and area coverage.Clearly, the noise may affect the ice dynamics, and a user must be aware of that.Fully eliminating it could be both difficult and expensive in terms of CPU time.Bouillon et al. (2013) proposed a modified EVP approach in which subcycling is fully detached from the physical time stepping.It can be considered as a pseudo-time solver for the VP rheology.In this case one writes for stresses and for the velocity.Here α and β are some large constants.The superscript p is related to pseudo-time iterations, replacing the subcycling of the standard EVP, and n is the index of external time stepping.Fields are initialized with values at time step n for p = 1, and their values for the last iteration p = N EVP are taken as solutions for time step n + 1.In order that CFL limitations are satisfied, the product αβ should be sufficiently large compared to π 2 P 0 t ( + min ) −1 m −1 x −2 (see Bouillon et al., 2013, andfurther comments by Kimmritz et al., 2015).The regime of the standard EVP scheme (N EVP = 120 and T = t/3) will be approximately recovered for α = β = 80 (for σ 1 ) and N EVP = 120, but much larger values have to be used on fine meshes to warrant the absence of noise in strain rates and viscosities.The stability requirements here are very similar to those of the standard EVP method if expressed in terms of N EVP and, likewise, become more restrictive for finer meshes.For numerical convergence, N EVP should exceed α and β (for the same reason that T is a fraction of t in the standard EVP).
One expects that if this scheme is stable and converged, it would produce solutions identical to those of a converging VP solver, while the standard EVP scheme may slightly deviate.We will return to this in Sects.3.4 and 3.5 where the time stepping is discussed.In practice, it will seldom be run for full convergence, which is rather expensive, and some difference will be preserved.
FESIM implements the three approaches mentioned above, which will be referred to further as VP, EVP and mEVP.The reason for keeping all of them is twofold.First, it facilitates the comparison of results with other models which may use one of these approaches.Second, their numerical efficiency and performance depend on applications, and one may wish to select the most appropriate one for a particular application.

Numerical implementation
We first describe spatial discretization, and then the discretization in time.

Finite-element discretization of ice transport equations
This section explains the FE spatial discretization, which is based on linear continuous functions defined on triangles.
The original motivation for this choice was the ability to share the infrastructure with the ocean model, which is based on the same discretization.The transport Eq. ( 7) are solved in two steps: first, scalar quantities are advanced with the right hand sides (rhs) of tracer equations set to zero.Then tracers are updated with account for thermodynamic sources and sinks in a pointwise manner.We therefore limit ourselves to homogeneous equations.In what follows, the superscript n will denote external time steps and p subcycling time steps in solvers, as in the discussion above.Subscripts j and k will denote quantities related to vertices (nodes) of triangular mesh.It is hoped that they will not be mixed with the notation for coordinate indices of tensors.For the mesh indices the agreement on summation over repeating indices will only be kept for matrix-vector products.
The tracer equations are solved with the FE Taylor-Galerkin (TG) method (see, e.g., Zienkiewicz and Taylor, 2000, p. 47), which is analogous to that of Lax-Wendroff for FV.One writes for the concentration and substitutes and In the last case the velocity is considered steady during the tracer time step.This still provides the second order in time if velocity and tracers are considered to be shifted by a half time step (asynchronous time stepping).The resulting equation provides the second order in both time and space (for linear functions).Here G is the modified flux vector, with a diffusive flux that exactly compensates for the first-order error in the time derivative.Note that it does not introduce dissipation.The ice and snow thickness equations are solved similarly.
To solve the tracer Eq. ( 20) with the FE method one first projects it on an appropriate set of test functions M j , and then integrates it by parts to obtain where is the boundary of the domain S. At the solid boundary (G • n = 0) or an open boundary located far from the icecovered region (so that a = 0), the boundary integral is zero.We will assume that this is the case.The procedure outlined above gives the equation in a socalled weak form.The discretization is obtained by expanding scalar fields and velocities into series: and similarly for h ice , h s , and components u and v of the velocity vector u.We use continuous Galerkin discretization implying that M j = N j , and that functions N j are continuous across the boundaries of triangles.We select N j as a linear function associated with vertex j of the triangular mesh.
It equals one at vertex j and decays linearly to zero at all neighboring vertices; the expansion above is simply the linear interpolation and summation is over all vertices.As a result, the Galerkin system of equations on nodal values of ice concentration a k (same for (h ice ) k and (h s ) k ) is obtained: where Note that summation is implied over k (matrix-vector product).It will be reminded in some cases below too.A similar procedure is used to obtain discretized momentum equations.The mass matrix M j k is not diagonal, but has a limited bandwidth (defined by the number of neighbors).Its appearance is what makes the method different from the FV Lax-Wendroff implementation.Indeed, it is easy to check that the latter would lead to the same result on mediandual control volumes (obtained by connecting triangle centroids with mid-edge points), but with the diagonal lumped mass matrix M L j k , whose diagonal entries are sums of rows of M j k , and other entries are zeros.Two points should be mentioned here on practical implementation.First, the velocity field is linear on triangles, so computations of operator A j k should be formally done with account for this.Doing so would not, however, improve accuracy compared to just using mean velocities on triangles, which simplifies computations.Second, true iterative solution of equations involving mass matrices, written schematically as M j k b k = c j , is expensive and is never attempted.Instead, one does three it- Doing more iterations does not improve dispersive properties of the method, yet doing just one (lumping) deteriorates the method rather noticeably.
The presence of a consistent mass matrix in the TG method effectively removes a significant portion of dispersion related to the Lax-Wendroff method.However, remaining dispersive errors may still be damaging.For this reason, the approach is augmented to the FE-FCT method as proposed by Löhner et al. (1987).In this method, the TG solution above serves as the high-order one, and will be denoted as a n+1 k .The loworder solution a n+1 k is obtained by adding artificial dissipation to the rhs and replacing the consistent mass matrix with the lumped one on the left hand side (lhs), which leads to a monotonic solution provided the parameter γ FCT is sufficiently high (about 1).The difference between the high-order solution a n+1 k and the monotonic low-order solution a n+1 k is due to the antidiffusive flux contribution: The rhs of the last expression is split into contributions from separate elements.They are limited as detailed in Löhner et al. (1987) and assembled back to recover a monotonic solution a n+1 k instead of a n+1 k .By construction, the solution method is conserving.Indeed, because j N j (x, y) = 1, j A j k = 0, and j M j k a k is the area integral.Also, j M j k a k = j M L j k a k , so that the simple iterative procedure above preserves conservation.According to Budgell et al. (2007) the FCT method shows second-order convergence in simple advection tests.Note, however, that the ice velocity is divergent and may thus lead to the formation of local extrema in scalar fields.The FCT scheme may therefore result in excessive smoothing of extrema.Yet it does so for the antidiffusive fluxes only, the low-order solution will react to the divergence of the velocity field.For this reason the parameter γ FCT should be taken at minimum compatible with stability and preservation of positivity.
Despite the fact that the FCT limiting doubles the computational cost of advection (compared to using solely the TG method), the burden remains small compared to the cost of solving for ice velocities.

Computation of strain rates and stresses
Similar to the thicknesses and concentration, ice velocities are considered to be linear functions on elements: The strain rates are therefore elementwise constant.At this point we need to take into account sphericity and peculiarities coming from the derivatives of metric terms.We use the spherical coordinate system with poles at land to avoid the pole singularity.In spherical coordinates (φ, θ ) and Here R it the Earth radius.We approximate the geometry as locally flat on triangles, which is equivalent to replacing cos θ in (1/ cos θ )∂/∂φ by its estimate on elements.If we use a local Cartesian frame of reference on each element with the x and y axes oriented along the directions e φ and e θ , we can then write ∂ x and ∂ y instead of (1/R cos θ )∂/∂φ and (1/R)∂/∂θ , respectively.With the same accuracy we make an elementwise-constant estimate of the metric differentiation term, so that the expressions above become where m f = tan θ/R is the metric factor.These expressions for the strain rates are further used to compute the components of stresses which would then be naturally treated as elementwise constant too.Although the ice strength P would be more naturally modeled as a linear function because h ice and a are represented in that way, the estimate of the ice strength gradient at vertex points will be the same if P is averaged to triangles, i.e., treated as elementwise constant.To further simplify computations we estimate h ice and a on triangles as the mean over vertices.This makes all components of stresses elementwise quantities, so that time stepping of stresses in EVP and mEVP becomes an algebraic operation on triangles.Formally projecting the last equations on functions M c = 1 on triangle (cell) c gives and Here summation is over vertices k of cell c, hence the symbolic notation k(c).The expression for the ice strength is computed as . With the strain rates and ice strength known, and the stress components are easily computed on elements.

Spatial discretization of momentum equation
Rigorous finite-element implementation of the momentum equation would involve mass matrices and would be too time consuming in the case of EVP and mEVP solvers.For that reason some simplifications are required.Luckily, mass matrices are not important here, as no compensation of discrete errors can be achieved with their help.We therefore use nodal quadratures in all terms that do not involve spatial derivatives.Multiplying Eq. ( 1) with test functions, integrating over the domain, and integrating the rheology term by parts, one gets Here N j is a shortcut for either (N j , 0) or (0, N j ), so that Eq. ( 34) is a set of two equations obtained by projecting on x and y directions, the second term on the rhs involves the dyadic product of two tensors and the last term involves the contraction of the stress tensor with the unit vector normal to the boundary.On substituting the expansions in N k for velocities, we approximate the lhs of Eq. ( 34) as where m k = ρ ice (h ice ) k + ρ s (h s ) k and M L j k is a shortcut for two "vectors" (M L j k , 0) and (0, M L j k ).Similarly, the first term on the rhs is Summation over k implied in these equations is trivial because the lumped mass matrix is diagonal.The entries of the diagonal lumped mass matrix (for j = k) are just the areas of median-dual control volumes associated with vertices, i.e., one-third of the sums of areas of triangles containing the vertex considered.
The second term on the rhs of Eq. ( 34) leads to the following contributions to equations for local x and y directions: Here c(j ) are the indices of cells containing vertex j (spanned by test function N j ) and A c is the area of cell c.Notice that, because of metric differentiation, applying ∇ to any of (N j , 0) or (0, N j ) also gives a contribution projecting on the other vector.
In the third term on the rhs of Eq. ( 34) computations of the slope term are simpler because the gradient of scalar field H does not involve differentiation of metrics.We use the nodal quadrature for the mass, which results in with summation over k implied.Here G x j k = N j ∂ x N k dS and similarly for the y-equation component.Clearly, and likewise for the y equation.
The last term in Eq. ( 34) involves only vertices j on the boundary.We do not need equations there in the no-slip case, which is used by us, because zero velocity will be prescribed by the virtue of boundary conditions.Leaving equations there but omitting the tangent component of this term would impose free-slip boundary conditions.

Time stepping and the implementation details of VP solver
As mentioned above, large values for viscosities in the VP case would lead to severe CFL limitations in the case of explicit time stepping.This suggests to account for the stress term in the ice motion equation implicitly: However, since the viscosities in σ are functions of the velocity field, the expression for σ should be linearized (by estimating viscosities on time step n) in order to use standard iterative solvers.The "implicitness" is recovered by doing (Picard) iterations, when the velocity of the previous iteration is used to estimate the viscosities for the current iteration.Note that friction between ice and ocean is linearized and taken implicitly too.
This approach is suboptimal because of the need to solve a problem for a matrix of dimension 2N , where N is the number of surface nodes (vertices).The nonzero entries in each row come from both u and v contributions in this case, which would make matrix-vector multiplications more expensive too.
The now traditional way of handling this problem was proposed by Zhang and Hibler (1997).In that case one makes implicit the terms involving u in the x equation and terms involving v in the y equation.This still requires assembling two matrices and preconditioning them.The approach employed by us was formulated by Hutchings et al. (2004).It is similar in spirit to that of Zhang and Hibler (1997), but allows us to use the same matrix for u and v.This considerably reduces the computational cost if general-purpose iterative solvers (like PETSc) are used.Its essence is in writing the stress tensor (3) in the form and making implicit only the first term on the rhs of this expression.Since the eigenvalue of the implicit operator is larger in this case than in the algorithm of Zhang and Hibler (1997), the method is stable.Yet its convergence rate is not necessarily better because it introduces an artificial residual through ζ (∂ i u j ).The rest of the implementation resembles that of Zhang and Hibler (1997).It consist of three steps.
The first two of them are iterations of the scheme where, as above, p is the index of iterations, and n of time stepping.In the original procedure p = 1, 2, but (Picard) iterations can be repeated to arbitrary high p = N p .For p = 1 the superscript * implies that the quantity is estimated at time step n.For p = 2, u * = (u p−1 + u n )/2, F * = F (u * ), and same for the viscosities on the lhs, following Zhang and Hibler (1997).For p > 2 (if N p > 2) the starred quantities are those at iteration p − 1.In the expressions above, F denotes the explicit part of the ice reaction.The final (third) step updates the Coriolis term to the implicit one: Because of the need to keep the same matrix in u and v equations, the terms associated with metric differentiation in the lhs operator are all put on the rhs (added to those of F ), and the discretization of the operator part is straightforward.For convenience, we write down F in the finite-element discretization.We first omit the terms arising from metrics differentiation, as they are more conveniently taken into account separately all together.Since and the divergence of the stress tensor multiplied with test function N j and integrated by parts will lead to the following contributions to the rhs of the u and v components of the momentum equations: All derivatives and P are elementwise constant, so the integrals are equivalent to summation over the cells spanned by N j .
It is easy to see that all "metric differentiation terms" lead to the additional contributions and respectively, to u and v equations.The last terms in both contributions require integration of test functions, which gives A c /3 on each cell involved.The operator matrix is assembled in the standard sparse format on each time step.In order to reduce the computational load in the course of iterative solution, the matrix entries in the rows corresponding to nodes where the ice concentration is less than a small critical value are set to one at the diagonal, and zero otherwise.The rhs vector is corrected accordingly, and set to zero (default) or to the ocean velocity or to the velocity of the previous time step.The PETSc solver with ILU (incomplete lower-upper) preconditioning is used to solve the resulting matrix problem.
In theory, the tolerance does not necessarily need to be very small as the solution procedure is repeated on every time step, and the solution cannot diverge very much from the previous solution.However, on unstructured meshes a small tolerance can sometimes be required to achieve an acceptable accuracy on elements of differing size.Also, higher solver accuracy can be needed in quasistationary regimes, to properly handle areas where is small.Our experience with PETSc is that while a tolerance of 10 −6 may be sufficient on relatively uniform meshes, it should be at least two orders of magnitude smaller if the size of mesh elements varies by a factor of 5 or more (see also discussion of convergence below).
There is always some sensitivity to the mesh, domain geometry and preconditioning; users are advised to experiment with the available options of the solver.
The expressions for the two last terms have been given above (Eqs.37, 38) and M j = M L jj with no summation (it is the area associated with vertex j ).The fields are initialized with values at time step n.Pseudo-time stepping of the momentum part of mEVP is given by Eq. ( 16) with the terms interpreted similarly as in the equations above.Now, when all equations are written, we can discuss the differences between the methods.The differences between the EVP and mEVP are subtle (apart from the difference in variables used to organize subcycling).First, (i) as can be seen comparing Eqs. ( 10)-( 12) with Eqs. ( 13)-( 15), the EVP uses different rates for σ 1 on one hand and σ 2 with σ 12 on the other to approach the VP rheology.Second, (ii) after N EVP iterations are done, the EVP scheme estimates the time derivative of velocity based on the last substep, while mEVP employs the estimate over the entire time step t.Third, (iii) there is damping in mEVP introduced by β, which helps to equilibrate the solution over the places where ice is weak.
One does not expect large discrepancies between both approaches.However, it turns out that (i) has a negative impact on stability (cf.Bouillon et al., 2013), which is why mEVP is more robust, as will be demonstrated below.At the end of the external time step the VP and mEVP solutions satisfy the same equations.To summarize, all three methods are expected to behave approximately similar, and the main point is the convergence of their solutions (and hence stability).

Box test case
The model described above is routinely used with FESOM both in an ice/ocean-only version or in a version coupled to an atmosphere model, so that its practical performance can be judged by the results of respective papers (see, e.g., Sidorenko et al., 2011;Wang et al., 2012;Wekerle et al., 2013;Timmermann and Hellmer, 2013;Wang et al., 2014;Sidorenko et al., 2015) and is not repeated here.Thus far FE-SOM was run only with the EVP solver (since it was the first one available) and the comparison of the performance of the three available versions in the global setup is the subject of future work.Here we will use a box test case without thermodynamic forcing, with an intention to demonstrate similarities and disparities in the performance of VP, mEVP and EVP algorithms.This will be more difficult for realistic simulations where many other factors may contribute.
The setup follows that used by Hunke (2001), with the difference that islands are removed, geometry is spherical and the mesh is an unstructured one with variable resolution as used in Losch and Danilov (2012).The square box is of approximately 11 • by 11 • in size (with the side lengths L x and L y ) and the resolution is varied approximately from 40 to 10 km from the south to the north, as shown in Fig. 1.It will be seen below that noise, if excited, appears at the fine mesh part, as could be anticipated.Apart from this, no other implications of mesh unstructuredness will be mentioned here to keep the discussion concise and concentrated on the algorithm performance issues.
Ice is driven by the wind stress τ = C a ρ a u a |u a |, with C a = 0.00225.Here ρ a is the air density and the wind velocity (in m s −1 ) is taken as u a = 5 + (sin(2π t/T ) − 3) sin(2π x/L x ) sin(πy/L y ) and v a = 5 + (sin(2π t/T ) − 3) sin(2πy/L y ) sin(π x/L x ), where T =4 days.The ocean velocity (in m s −1 ) is selected as and the elevation H is computed by geostrophy.The coordinates x, y are the longitude and latitude counted from the southwest corner of the box.The ice thickness is 2 m initially and the ice concentration grows linearly from 0 to 1 in the west-east direction.The results of simulations at the end of 1 month are shown.
We start by comparing VP and mEVP solutions.In case A advection is switched off, and we compare the convergence of solutions obtained with different methods.In cases B and C the advection is switched on, they differ by the value of min : 2 × 10 −9 s −1 (B) and 2 × 10 −11 s −1 (C).Case A takes min of case B. Figures 2 and 3 show, respectively, the zonal velocity and (upper left panels) and the differences in solutions obtained by different methods in case A. We take the mEVP solution with α = β = 500 and N EVP = 1000 as a reference one (mEVP500), for modifications seen in runs with higher α, β and N EVP are very small.The other solutions shown are those obtained with mEVP, but α = β = 250 and N EVP = 250 (mEVP250), and with VP, but in the regime with 2 (VP2p) and 10 (VP10p) additional Picard iterations (which means that N p = 4 and 12, respectively).It is immediately seen from the velocity comparison that mEVP250 is far from convergence (there is a large-scale pattern in the velocity difference) and that it contains noise in the field of .Note that the noise is seen over the fine part of the mesh, as stressed in Losch and Danilov (2012), because it is more difficult to satisfy the stability requirement when the mesh is refined.So the parameters of the mEVP and the number of subcycles should be sufficiently large.Note that the same is also true for the standard EVP.The traditional practice of running it with relatively low subcycling numbers (N EVP = 120 is commonly used) may lead to noise in over places where it is sufficiently small.
The difference between the two VP solutions and mEVP500 is much smaller and is largely concentrated at the front between the moving and nearly stopped ice.However, one sees that there is a basin-scale pattern in the velocity difference in the bottom left panel of Fig. 2, which is the indication of the lack of convergence of the VP solution over the area where ice is moving.Indeed, it almost disappears on increasing the number of Picard iterations (bottom right panel).Simultaneously we see the substantially improved agreement Ice zonal velocity (m s −1 ) in case A (advection is switched off) after 1 month of simulations in mEVP500 (top left) and differences between the solutions obtained by different methods: mEVP250-mEVP500 (top right), VP2p-mEVP500 (bottom left) and VP10p-mEVP500 (bottom right).mEVP250 does not converge, and VP2p is closer to convergence but still with noticeable errors.Additional Picard iterations in PV10p substantially reduced differences between the mEVP and VP solutions.
between the patterns of in Fig. 3.The remaining discrepancy is due to errors in both, EVP500 and VP10p, solutions, eliminating it will require increasing the number of subcycling steps and iterations even further and is not pursued.We conclude that mEVP and VP converge to each other if one takes care that both are sufficiently accurate.
Reaching full agreement between mEVP and VP solutions is more difficult if the ice advection is on, because errors may accumulate in this case with time.Smaller values of min additionally complicate the issue.In the presence of advection, ice is pressed into the northeast corner of the mesh, piling up there.The western part of the basin becomes free of ice, so that there are two fronts no ice-moving ice and moving ice-nearly stopped ice.We concentrate on the differences in the northeast corner, errors along the fronts depend on minor details and are difficult to eliminate.
The results of case B are given in Figs. 4 and 5 which present h ice and , respectively, after 1 month of model time.
Here we compare three VP solutions with the mEVP500 reference simulation.We checked that increasing α and β to 1000 with subsequent increase of N EVP to 2000 in mEVP does only small changes to the field of compared to those seen for the VP solutions.The solution labeled VPb is obtained with the basic algorithm (N p = 2), and VP10p and VP25p correspond to using 10 and 25 additional Picard iterations, respectively.While the difference in ice thicknesses remains small and is only slightly affected by the number of iterations in the VP solutions (patchiness in the difference  panels of Fig. 4 is due to the finite accuracy of output), there is substantial improvement in the correspondence between the mEVP and VP solutions for as the number of Picard iterations is increased.The fact that the differences in the ice thickness in the northeast corner stagnate hints that they evolved from some minor implementation details.Since the total ice volume is conserved, these errors are connected to those in the front position.They are rather small to be of practical importance.(s −1 ) after 1 month of simulations in case B in mEVP500 (top left) and differences between the solutions: VPb-mEVP500 (top right), VP10p-mEVP500 (bottom left) and VP25p-mEVP500 (bottom right).Additional Picard iterations in the VP method lead to substantially reduced differences between the solutions in the northeast corner.VPb reproduces a much stronger (smaller ), but additional Picard iterations make it weaker and closer to mEVP500.Finally,case C (Figs. 6,7) shows that reaching agreement between the mEVP and VP runs for a much smaller min requires an even larger number of Picard iterations (and also more subcycling in mEVP, although the improvements seen are less substantial).The mEVP500 solution in this case contains some noise in , and is replaced by mEVP1000 obtained with α = β = 1000 and N EVP = 2000.We also consider the standard VPb solution and the solutions obtained with 100 (VP100p) and 200 (VP200p) Picard iterations.As in case B, the Picard iterations do not change the difference between ice volumes very much, but have substantial impact on the field of .Similarly, VPb produces a stronger ice (smaller ) in the northwest corner, which is partly made weaker by increased number of Picard iterations.Of particular interest is the structure in the compression zone of VP solutions, which is sensitive to the number of iterations.There is some sensitivity of band structure to the change of solver tolerance and time step.This hints that one deals here either with incomplete convergence or some internal instabilities in the iterative procedure, a question we postpone for the future.We see that it is much more difficult to minimize the difference between mEVP and VP solutions if min is taken smaller.
Since the intention of min is to provide regularization, it should not be made excessively small unless there is motivation for that.
The next pair of figures (Figs. 8,9) compares the performance of EVP and mEVP solvers.We use min = 2 × 10 −9 s −1 , and three EVP solutions denoted EVP3_100 ( t/T = 3, N EVP = 100), EVP3_500 (N EVP = 500) and a special solution, EVP4_1000 ( t/T = 4, N EVP = 1000), obtained by removing e 2 from the second terms on the left hand side of Eqs. ( 11) and ( 12) and putting it to the denominator of the right hand side.After this manipulation EVP becomes almost identical to mEVP (all components of the stress tensor decay to the VP limit at the same rate), except for the differences in the velocity time stepping.In this case one may identify α with 2T N EVP / t.Solution EVP3_100, with parameters typical for large-scale applications, shows noisy over the area with compressed ice.Increasing the number of subcycle steps improves the agreement (Fig. 9, bottom left) but it still remains noisy.The noise takes the form of a wave structure.Simulations with further increased N EVP (1000 and 2000, not shown) improve the agreement, but only slightly.Similarly, varying T is of no avail.However, the situation improves dramatically if the decay rates for stresses in Eqs. ( 10)-( 12) are made similar, as indicated by the bottom right panel in Fig. 9.The noise disappears.While the remaining discrepancy in over the stiff ice can be further reduced, some differences will persist because of the different treatment of the momentum equation.The central circular spot over weak ice is common to all three solutions.Here the contribution from rheology is not dominant, and the difference is entirely due to the time stepping of the momentum equation.We therefore conclude that it is the difference in the damping rates in the equations for stresses, Eqs. ( 10)-( 12), in the standard EVP which is the main factor responsible for the noise seen in the field of .More detailed (s −1 ) after 1 month of simulations in case C in mEVP1000 (top left) and differences between the solutions: VPb-mEVP500 (top right), VP100p-mEVP1000 (bottom left) and VP200p-mEVP1000 (bottom right).Additional Picard iterations in the VP method substantially modify the differences, reducing them in the northeast corner.The convergence is not reached even for 200 Picard iterations.
analysis of this statement is needed.If we now turn to the patterns of ice thickness, we see that even in EVP3_100 and EVP3_500 with noisy the simulated mean ice thickness agrees rather well with the mEVP solution, with differences of about 10 cm at maximum.The difference virtually disappears for the special case of EVP4_1000.
In summary, given the sensitivity of the field of to the solution procedure, one should be cautious to discuss its detail unless the convergence has been tested.Judged from this perspective, the VP and mEVP approaches provide more consistent behavior than the EVP.However, even with them, one should realize that there might be some sensitivity to the implementation detail.For example, the VP solutions discussed here have been obtained with a tolerance of 10 −8 in the PETSc solver; using a tolerance of 10 −6 leads to changes in comparable in magnitude to the effect of varying the number of Picard iterations.We have not seen benefits from making the tolerance even smaller (10 −10 ), but this may change in other applications.Additionally, there is some sensitivity to the time step interval t.Finally, the lack of full agreement in the pattern of in VP and mEVP simulations, especially for the low min = 2 × 10 −11 s −1 in case C, can partly be due to the particular implicit/explicit splitting of the stresses, and we cannot exclude that the original splitting of Zhang and Hibler (1997) will converge somewhat differently.Note that the mEVP method shows less sensitivity to details than the VP method if α and β are sufficiently large to ensure the absence of noise in the solutions and if N EVP is sufficient for convergence.The ice mean thickness and concentration, in contrast, show a much more robust behavior, and are much more consistent, even in the presence of noise in .Still, the presence of noise pushes simulations on a dangerous ground and should be avoided.In many practical cases the VP, mEVP or EVP solvers will be run in a "partially converging mode" when accuracy is achieved over a number of steps under conditions that forcing does not change much over a time step.Numerical stability and lack of noise (for the EVP and mEVP methods) will remain an issue to pay attention to.

Numerical aspects: spatial discretization
The finite-element discretization of sea ice dynamics employed by FESIM works in a robust way on unstructured triangular meshes.We now discuss how it relates to other unstructured-mesh discretizations proposed in the literature.
We first note that the FE P 1 − P 1 implementation can easily be cast in a FV form as explained in the Appendix.As concerns the purely dynamical (momentum) part, there is almost no difference in the final result to the FE discretization because of the lumping of the mass matrices we use for dynamics.One may wish to select a transport scheme that differs from FE-FCT, but the only motivation behind this can be the availability of a more accurate and efficient FV scheme.Our experience is that reaching the accuracy of the FE-FCT scheme would require a better than third-order method in the respective FV-FCT algorithm.As mentioned above, the presence of a consistent mass matrix in the FE transport equation efficiently compensates for a significant part of dispersion, which explains its good performance.
The vertex placement of variables we used is an analogue of the A grid in the traditional (Arakawa) terminology.A different A-grid implementation with the cell (triangle centroid) placement of variables was proposed by Hutchings et al. (2004).The discretization is straightforward if the FV approach is used and if the velocity derivatives on each triangle are computed by, e.g., the least square fit using the velocities on this and three neighboring triangles.The potential problem of the cell-based placement is a somewhat unfavorable stencil used in the computation of stress divergence.Indeed, it involves not only the nearest neighbors, but the neighbors of neighbors.We therefore consider the vertex placement of variables to be an easier choice.
The implementation adopted by FVCOM (Gao et al., 2011) is also a FV one, with velocities placed at cells and scalars at vertices.We tested this variable placement while developing the sea ice model to complement the ocean circulation model based on the staggered cell-vertex discretization.Because of an excessively large velocity space (the number of triangles is approximately twice that of vertices) it is prone to noise in velocities along the ice edge and was therefore abandoned in favor of the vertex-vertex scheme.Once again the vertex placement of velocities and scalars seems to be a more robust option.
Finally, the discretization proposed by Lietaer et al. ( 2008) is a FE one, based on nonconforming linear functions to represent the velocity vectors, with velocity degrees of freedom placed at the edges and elementwise-constant scalars.It also has a too-large velocity space and is not optimal in this re-spect.Additionally, placing scalars at centers would be suboptimal for representing the ice strength gradients: a nonconforming linear function used for velocity spans only two elements with a common edge, and two ice strength values at centroids give only one component of the gradient.
Thus, despite its simplicity the discretization in FESIM deserves attention as a balanced choice.Work is planned on augmenting it with a multi-category ice functionality.

Numerical aspects: VP/EVP convergence
There is ongoing discussion on the convergence of traditional implementations of VP and EVP, with indications that convergence is lacking (see, e.g., Lemieux and Tremblay, 2009;Lemieux et al., 2012).It partly motivated the development of new approaches such as the Jacobian-free Newton-Krylov solver (see Lemieux et al., 2012), which intends to improve the convergence of the VP method, but is too CPUdemanding, and also served as a motivation behind the new formulation of EVP in Bouillon et al. (2013), referred to as mEVP here.However, Bouillon et al. (2013) mention that they fail to reach converging mEVP solutions.The analysis of Kimmritz et al. (2015) shows that mEVP does provide converging solutions, but only when α and β are sufficiently large and N EVP is larger than any of them.From the theoretical viewpoint the mEVP and VP methods should lead to identical solutions if converged, and the solutions obtained with EVP may slightly deviate from them.The box test cases above illustrate that the solutions can be made rather close, but reaching full agreement between them might be too expensive computationally and require adjusting minor details of the algorithms.
The stability (and convergence as a result) of (m)EVP solvers is sensitive to the mesh size, and will generally deteriorate if the mesh is refined.Larger α, β, N EVP are to be expected on finer meshes, and it is the user's responsibility to select values providing the absence of noise in the fields of divergence and .Note that the issues mentioned here are in full measure relevant for other models, including those formulated on structured meshes.While in realistic applications they can be hidden behind much larger uncertainties in parameterizations of mechanical and thermodynamical forcing, one should be sure that the dynamical operators the model relies on behave in a predictable and understandable way.

Practical aspects: CPU load
Computations of stresses and their contributions to the rhs of momentum equation are rather expensive in models formulated on unstructured meshes (compared to their structuredmesh counterparts) mainly because of the lack of directional splitting and, in the case presented, also because the number of triangles is twice as large as the number of scalar degrees of freedom.For this reason, one computation of the rhs (done N EVP times per external time step in EVP and mEVP solvers) is substantially more expensive than one matrix-vector multiplication in the iterative matrix solver in the VP method.On the other hand, the number of iterations needed to reach convergence to the specified tolerance may depend on the ice distribution and domain geometry and the number of required Picard iterations can be high.One has to take into account the time spent on assembling the stiffness matrix and preconditioning it.Any comparison is even more complicated because full convergence of mEVP and VP methods will not necessarily be attempted in practice.For this reason no general recommendation can be given here.Just for information, we present the results for the box test case above: the time step of mEVP500 with N EVP =1000 takes 0.55 s on eight cores of an old IBM BladeCenter JS22, to be compared with 0.88 s for VP25p and only 0.065 s for VPb, and there is approximately linear dependence on N EVP and the number of Picard iterations N p .Since VPb (N p = 2) provides a very reasonable solution for the ice mean thickness and since the field of , despite the lack of convergence, is smooth in this case, it can still be used and will be a faster option than mEVP500 with N EVP = 1000.They will be close to each other if we run mEVP500 with N EVP = 120, sacrificing convergence but keeping stability.As mentioned, the comparison in a realistic global configuration is the subject of future work.

Conclusions
FESIM, the sea ice component of FESOM v.1.4,is described here.We focus on the dynamical part of the model in this documentation.The new EVP solver (mEVP) proposed by Bouillon et al. (2013) leads to solutions approaching those of the VP solver if both are run toward convergence.However, it is expected that some differences between their results would still persist in practical usage.While the mEVP Eqs. ( 13)-( 16) algorithm shows better stability in our tests than the standard EVP algorithm Eqs. ( 49)-( 55), the performance of mEVP and VP is rather similar, and the CPU efficiency becomes the criterion to select between them.The mEVP method can be more convenient on massive parallel computers.As concerns the unstructured character of meshes, the implementation based on linear continuous elements is perhaps the easiest among the other possible choices.It shows robust behavior and serves well the tasks of multi-resolution modeling, as indicated by a growing list of practical applications using FESOM.An important issue to be kept in mind with respect to multi-resolution simulations is the sensitivity of stability and hence convergence to the mesh resolution.

Figure 1 .
Figure 1.Triangular mesh used in simulations.The resolution varies from approximately 40 to 10 km.Stability of EVP and mEVP on the fine mesh requires that α, β and N EVP be sufficiently large.

Figure 2 .
Figure2.Ice zonal velocity (m s −1 ) in case A (advection is switched off) after 1 month of simulations in mEVP500 (top left) and differences between the solutions obtained by different methods: mEVP250-mEVP500 (top right), VP2p-mEVP500 (bottom left) and VP10p-mEVP500 (bottom right).mEVP250 does not converge, and VP2p is closer to convergence but still with noticeable errors.Additional Picard iterations in PV10p substantially reduced differences between the mEVP and VP solutions.

Figure 3 .
Figure 3. Same as in Fig. 2 but for the "divergence" (s −1 ) after 1 month of simulations.Additional Picard iterations in the PV method lead to very good agreement between mEVP and VP solutions.

Figure 4 .
Figure 4. Ice thickness h ice (m) after 1 month of simulations in case B in mEVP500 (top left) and the differences between solutions obtained by different methods: VPb-mEVP500 (top right), VP10p-mEVP500 (bottom left) and VP25p-mEVP500 (bottom right).Additional Picard iterations in the VP method only slightly affect the differences.
Figure 5.(s −1 ) after 1 month of simulations in case B in mEVP500 (top left) and differences between the solutions: VPb-mEVP500 (top right), VP10p-mEVP500 (bottom left) and VP25p-mEVP500 (bottom right).Additional Picard iterations in the VP method lead to substantially reduced differences between the solutions in the northeast corner.VPb reproduces a much stronger (smaller ), but additional Picard iterations make it weaker and closer to mEVP500.

Figure 6 .
Figure 6.Ice thickness (m) after 1 month of simulations in case C ( min = 2 × 10 −11 s −1 ) in mEVP1000 (top left) and differences between the solutions obtained by different methods: VPb-mEVP1000 (top right), VP100p-mEVP1000 (bottom left) and VP200p-mEVP1000 (bottom right).The differences are small and additional Picard iterations only slightly change them.
Figure 7.(s −1 ) after 1 month of simulations in case C in mEVP1000 (top left) and differences between the solutions: VPb-mEVP500 (top right), VP100p-mEVP1000 (bottom left) and VP200p-mEVP1000 (bottom right).Additional Picard iterations in the VP method substantially modify the differences, reducing them in the northeast corner.The convergence is not reached even for 200 Picard iterations.

Figure 8 .
Figure 8. Ice mean thickness h ice (m) after 1 month of simulations in mEVP500 (top left) and differences between it and EVP solutions: EVP3_100-mEVP500 (top right), EVP3_500-mEVP500 (bottom left) and EVP4_1000p-mEVP500 (bottom right).The last EVP solution (bottom right), obtained with modified equations for stresses, shows the results almost identical to mEVP (see the text for details).

Figure 9 .
Figure 9. Same as Fig. 8 but for (s −1 ).Only the special solution obtained with the same decay rates in equations for stresses (bottom right) compares well to the mEVP solution.