1 Introduction
We consider the viscousplastic (VP) sea model, that was introduced by Hibler in 1979 and which is still one of the most widely used sea ice rheologies as detailed by Stroeve et. al. (2014). The model includes strong nonlinearities such that solving the sea ice dynamics at high resolutions is extremely costly and good solvers are under active research. Mostly, solutions to the VP model are approximated by iterating an elasticviscousplastic (EVP) modification of the model that was introduced by Hunke and Dukowicz (1997) and that allows for explicit subcycling. Alternatively the VP model is tackled directly with simple Picard iterations as described by Hibler (1979) or solved with Newtonlike methods as described by Lemieux et. al. (2010) or Mehlmann and Richter (2017b). All approaches are not satisfactory as they are extremely expensive and often are not able to give an accurate solution in reasonable computational time. It is therefore of utmost important to reduce the complexity of the computations, e.g. by using coarse meshes and large time step sizes, as long as this does not deteriorate the accuracy assumptions.
We derive an error estimator that identifies the errors coming from spatial and temporal discretization. Furthermore, the error estimator allows for a localization of the error to each element and each time step such that local step sizes can be adjusted. This goal oriented error estimator for the viscousplastic sea model is an extension of the dual weighted residual method that was introduced by Becker and Rannacher (2001). The aim of the estimator is to identify discretization errors between the unknown exact solution and the numerical approximation , where indicates the temporal and the spatial discretization parameter, in functionals . These functionals can be any measures of interest, e.g. the average sea ice extent in a certain time span
(1) 
We denote by the ice concentration (one component of the solution which will be introduced later), by the spatial domain of interest and by the time span of interest, e.g. the summer months. The error estimator will give approximations to which can be attributed to a spatial error, a temporal error and to a splitting error  coming from partitioning the system into momentum equation and balance laws, which is the standard procedure in sea ice numerics, see Lemieux et. al. (2014). The estimation of errors in space and time for parabolic problems was discussed by Schmich and Rannacher (2012). Here, we extend the application to the VP model and additionally consider the splitting error. Lipscomb et. al. (2007) pointed out that decoupling the system in time can lead to a numerical unstable solution such that a small time step is required to achieve a stable approximation. Lemieux et. al. (2014) introduced an implicitexplicit time integration method, which resolves this issues and allows the use of larger time steps. The error estimator will be able to predict the accuracy implications of this temporal splitting.
The convergence of approximations of the sea ice velocity on different spatial resolutions is studied by Williams and Tremblay (2018) for a one dimensional test case. The authors observe that the simulated velocity field depends on the spatial resolution and found that the mean sea ice drift speed rises by 32% by increasing resolution from 40 km to 5 km. The temporal and spatial scaling properties of the mean deformation rate and the sea ice thickness are studied by Hutter et. al. (2018).
The paper is structured as follows. In Section 2 we start by presenting the sea ice model in strong and variational formulation which is required for the Galerkin finite element discretization in space and time. Further give details on the partitioned solution approach. In Section 3 we derive the goal oriented error estimator for the sea ice model and describe its numerical realization. We numerically analyse the error estimator in Section 4 and conclude in Section 5. For better readability we keep the mathematical formulation as simple as possible and refer to the literature for details.
2 Model Description and Discretization
Let be the spatial domain. We denote the time interval of interest by . The dynamics of sea ice is described by three variables, the sea ice concentration , the sea ice thickness and the sea ice velocity , such that the complete solution is given by . The VP sea ice model as introduced by Hibler (1979) consists of the momentum equation and the balance laws
(2)  
with and . The forcing term models ocean and atmospheric traction
with the ocean velocity and the wind velocity . By we denote the ice density, by the Coriolis parameter, by the radial (
direction) unit vector. Following
Coon (1980) we have replaced surface height effects by the approximation . In this paper we focus on the dynamical part of the sea ice model such that we neglect thermodynamic effects and set and .Parameter  Definition  Value 

sea ice density  
air density  
water density  
air drag coefficient  
water drag coefficient  
Coriolis parameter  
ice strength parameter  
ice concentration parameter 
The system of equations (2) is closed by Dirichlet conditions on the boundary of the domain and initial conditions , and for ice thickness, concentration and velocity at time .
Finally, we present the nonlinear viscousplastic rheology which relates the stress to the strain rate
where is the trace. The rheology is given by
(3) 
with the viscosities and , given by and
(4) 
is the threshold that describes the transition between the viscous and the plastic regime. The ice strength in (3) is modeled as
(5) 
with the constant . All problem parameters are collected in Table 1.
2.1 Variational formulation and discretization
The dual weighted residual estimator by Becker and Rannacher (2001)
relies on a variational formulation of the system of partial differential equations in space and time. Galerkin discretizations are based on approximating the differential equation by restricting the search space of solutions to a finite dimensional space. In contrast to difference methods, where the discretization is based on an approximation of the differential operators, the equation itself is not changed. While the finite element method is a well established Galerkin method for the spatial discretization, time Galerkin discretizations are not widely used. Before describing the time discretization in detail we want to clarify that the Galerkin approach has to be considered mostly as a mathematical concept. The resulting discretization scheme is nearly identical to the backward Euler method and this is how we actually solve the problem. Nevertheless, the reformulation into a Galerkin approach is essential to derive the error estimate.
To start with, we multiply the three equations (2) with test functions , and and integrate over the temporal interval and the spatial domain . To shorten the notation we introduce the inner products in space and in the spacetime domain
Then, the variational formulation of the system of sea ice equations is given by the relation
(6) 
for all test functions , and , where are suitable function spaces.^{1}^{1}1This variational formulation is also called weak formulation, since it requires only first derivatives of the solutions and furthermore, since these derivatives are only required to be integrable functions.
The finite element discretization of (6) would consist of replacing the function spaces by discrete subspaces of finite element functions, e.g. , where is the space of piecewise linear functions defined on a triangulation of the domain .
In time, we proceed in a similar way and replace the solution and the test function by functions that are piecewise constant on the time partitioning
The time step size is denoted by . For simplicity we assume that is constant for all the time steps. Finally, the velocity is found in the spacetime Galerkin space
On each subinterval the velocity
(7) 
is a constant function in time with values in the discrete finite element space . Similar constructions are done for the ice concentration and the ice thickness . Since these discrete spaces are discontinuous in time (but the real solution is continuous in time), they are no subspaces of . The method is therefore called discontinuous Galerkin time discretization, in short: dG(0), see Schmich and Vexler (2008).^{2}^{2}2The “0” in dG(0) stands for piecewise constants
Similar to the finite volume method, which can be seen as the spatial counterpart to the dG(0) discretization we must add numerical fluxes in time to guarantee continuity of the solution in the limit . Here, we add jump terms at each discrete time steps
(8) 
where we denote by the rightsided value of the (possibly) discontinuous function , by the value from the left and by the jump at time . The real solution is continuous and it holds and . If the jump is not zero, the terms in (8) can be regarded as penalty terms.
To shorten notation we combine with and with and we assume that these functions come from function spaces and . The exact notation of all function spaces is introduced in Mehlmann and Richter (2017a). Then, the variational formulation (6) plus the additional jump terms (8) can be written in an abstract notation by introducing the form which simply collects all the integrals from (6) and the jumps from (8)
(9) 
The analytical solution is described by this variational formulations. All jumps vanish such that the solution is continuous in time. The discrete solution is given by restricting (9) to the finite dimensional discrete space
(10) 
In time, these discrete functions are piecewise constant, as defined in (7). In space, the discretization is quickly described: We define a conforming finite element space . In our implementation we use the space of piecewise bilinear functions defined on a quadrilateral mesh of the domain . Danilov et. al. (2015) consider standard linear finite elements on triangular meshes.
The goal oriented error estimator relies heavily on the concept of Galerkin discretizations: the discrete solution is found by limiting (9) to the finite dimensional space while the variational form itself is unchanged. This is in contrast to finite difference schemes that are based on discrete approximations of the differential operators. The most important tool in the analysis of Galerkin discretizations is the Galerkin orthogonality: the discretization error is orthogonal on the Galerkin space, i.e.
(11) 
which directly follows from subtracting (10) from (9), since is a subspace and the real solution satisfies (9) also, if tested with discrete functions .
2.2 Equivalence to the backward Euler method
The Galerkin formulation (10) appears rather abstract. We will show that it coincides with the backward Euler method that is most often used for solving the sea ice dynamics. On each interval the solution is constant in time, as defined in (7). As it is has discontinuities at and there is no natural coupling to the adjacent intervals. The only coupling between intervals is introduced by the jump terms (8), where couples to . These jumps take the role of the time derivative. There is no coupling of to or the previous intervals , as shown in Figure 1. Since is constant on it holds that . Further, all temporal integrals can be exactly evaluated with the box rule such that the spacetime Galerkin formulation corresponds to the timestepping scheme for
(12)  
Division by the step size reveals the classical backward Euler scheme which is standard in sea ice dynamics as described by Lemieux et. al. (2014). The Galerkin formulation (10) has to be considered as a mere concept of notation. Using the Galerkin formulation be essential for derivation of the error estimator.
Remark 2.1
The transport equations for and are under the constraints and . Without proper modeling of the thermodynamics there is no natural effect to guarantee these constraints. We will therefore apply a simple modification of the discrete system to take care of . We introduce the right hand side
that only gets active if and that will then force below one.
2.3 Partitioned solution approach
The discrete formulation (10) naturally splits into timesteps as shown in (12). The three components velocity , ice concentration and ice thickness however are coupled. It is standard to apply a partitioned solution approach in every time step, either by first solving the momentum equation for the sea ice velocity followed by the balance laws, or vice versa, see Lemieux et. al. (2014). We start with the momentum equation and replace the dependency of the ice mass and the ice strength given in (5) by the ice concentration and thickness from the last time step. To realize this decoupling within a Galerkin method we introduce the projection operator that projects (or , respectively) on the interval onto the value . For discrete functions this corresponds to and .
In terms of the backward Euler scheme, which is used for simulation, this means replacing and in the momentum equation (terms tested with ) by the values and from the previous step. For the Galerkin formulation this calls for a slight modification of the variational formulation , namely the introduction of the projection operator in the momentum part
(13) 
Once again, we indicate this variational form only for the formulation of the error estimator, the solution itself is computed by the backward Euler scheme (12) by replacing in the momentum equations by the previous approximations .
For the following we define the splitting error that is introduced by partitioning the sea ice system in time as the difference between both variational formulations
(14) 
3 Goal oriented error estimation
In this section, we derive a goal oriented error estimator for partitioned solution approaches. The new error estimator will be based on concepts of the dual weighted residual (DWR) method introduced by Becker and Rannacher (2001). The DWR estimator can be easily applied to all problems given in a variational Galerkin formulation and it has been applied to various problems for error estimation in space and time such as fluid dynamics (Schmich and Rannacher (2012)) or fluidstructure interactions (Richter, 2017, Chapter 8). The concept can also be extended to the estimation of modeling errors, see Braack and Ern (2003), an approach that has similarities to the treatment of partitioned solution schemes considered within this work.
Before deriving the error estimator for decoupled systems, we briefly recapitulate the idea of the dual weighted residual method.
3.1 The dual weighted residual method
To illustrate the idea of the DWR estimator we consider the heat equation in a spacetime variational formulation with
(15)  
Discretization is by means of the dG(0) method in time (which corresponds to the implicit Euler scheme) and finite elements in space such that we have to deal with the two problems
with the Galerkin subspace .
We start by defining an error measure, e.g. the average temperature at final time described by the error functional
Key to error estimation with the dual weighted residual method is the introduction of the dual solution , which is given as solution to the dual problem
(16) 
and its discretization
Computing all integrals, the dual variational formulation reveals the implicit Euler method of a problem that runs backward in time.
(17)  
(18) 
Instead of having an “initial value” at time we solve an initial problem (17) for . With the dual solutions and the error in the functional is given by
which directly follows by inserting in (16). This is allowed since comes from a Galerkin subspace. The error estimator is then completed by using Galerkin orthogonality (11
) to introduce an interpolation
of the dual solution and using that(19) 
This relation is called error identity. It cannot be exactly evaluated in practical application since it involves the unknown dual solution . However, there are well established strategies for approximating the local interpolation error and by reconstruction techniques, see Richter and Wick (2015). These interpolation errors act as weights that focus the residuals on areas with high sensitivity to the error functional.
The standard feedbackapproach for running adaptive simulations based on the DWR method is as follows:
Algorithm 3.1
Let be the initial time mesh, the initial spatial mesh. Let be the resulting spacetime function space.

Solve the primal problem

Solve the dual problem

Approximate the weights ,

Evaluate the error identity (19)

Stop, if is sufficiently small

Otherwise use the error estimator to adaptively refine the spatial and temporal discretization and restart with a finer function space on the refined mesh .
In Step 3, the approximation of the weights is the delicate part of the error estimator. Here, we must replace the unknown exact solutions and
. For heuristic approaches we refer to
Becker and Rannacher (2001), or in particular Schmich and Vexler (2008) or Meidner and Richter (2014) for time discretizations.Application of the DWR method will always require some numerical overhead, mainly by the computation of the auxiliary dual problem. It turns out that the dual problem is always a linear problem, also in the case of the fully nonlinear sea ice problem. In the following sections we describe the steps that are required for applying the DWR estimator to the sea ice model in a partitioned solution framework. We refer to Mehlmann (2019) for the full derivation.
3.1.1 Extension to the nonlinear case
For linear problems, the dual problem is just the transposed of the primal problem, compare (16). The general nonlinear setting requires a more involved approach. Here we cast the error estimation problem into the framework of optimization problems. We minimize the functional value under the constraint that satisfies the sea ice problem, . We introduce the Lagrangian
(20) 
where takes the role of the Lagrange multiplier. We then obtain the nonlinear error identity (since and )
(21) 
The right hand side is written as an integral
(22) 
and this integral is approximated with the trapezoidal rule
(23) 
with and . The third order remainder is omitted in practical application. The derivatives appearing in (22) are the directional derivatives and if we consider the definition of the Lagrangian (20) they take the form
Details are given in Section 3.2. To proceed with the nonlinear error identity (23) we now define and as the solutions to the linearized dual problems
(24)  
The specific form for the sea ice problem will be given in Section 3.1.2. Galerkin orthogonality (11) also holds for this linearized dual problem. Hence, the first term in (23) is zero and the complete nonlinear error estimator is given by
(25) 
where we used primal and dual Galerkin orthogonality (11) to replace the errors and by local interpolation errors. The first part is the primal residual and it is multiplied with the dual weights , the second part is the dual residual which is multiplied with the primal weights.
For details and many examples on the DWR method we refer to Becker and Rannacher (2001) or Schmich and Vexler (2008) who focus on the time dependent part. Algorithm 3.1 must not be changed for the general nonlinear case. Only the dual problem is constructed in a more complex setting. It is, however, still a linear problem.
3.1.2 The goal oriented error estimator for partitioned solution approaches
The partitioned solution approach is realized as a nonconsisting discretization of the sea ice equation: for discretization we replace the variational formulation by the splitting form . This means that the Galerkin orthogonality (11) is perturbed and it only holds
or, using the splitting error (14)
This requires subtle modifications in the derivation of the error estimate. Again, we only sketch the derivation and refer to Mehlmann (2019) for details.
Most important, the nonlinear error identity (21) must be based on two separate Lagrangians
where
(26) 
is based on the projection operator . Then, the DWR method is derived by dividing the error identity into the Galerkin error and into a splitting error by introducing
While the Galerkin part can be estimated as outlined in Section 3.1.1, the splitting part is quickly derived by using the definitions of the Lagrangians in (20), (26) and (14)
This error contribution can be evaluated since it only depends on quantities that are available, namely the primal and dual discrete solutions. We summarize:
Theorem 3.2 (DWR estimator for partitioned solution schemes)
Let and be primal and dual solutions to
Then, it holds that
(27) 
where is given in (23) and with the primal and dual splitting errors
By we denote an interpolation to the spacetime domain, by and we denote semidiscrete solutions which are discretized in time only.
Remark 3.3 (Weights)
The error estimator (27) depends on the unknown solutions and but also on semidiscrete solutions and , which are still continuous in space. All these terms must be approximated by suitable reconstruction techniques in order to evaluate the error estimator. In general, reconstruction is be a postprocessing mechanism: In time, we combine two intervals and reconstruct a linear function by connecting in with at , i.e.
(28) 
Then, we approximate . In space a similar procedure is done by combining the piecewise linear function on four adjacent quadrilaterals to one quadratic function. We refer to Richter and Wick (2015) and Mehlmann (2019) for details. The notation means: interpolation to linear, , functions on a mesh with double, spacing in time. Correspondingly, stands for the interpolation to the space of quadratic, , functions on the space with double spatial mesh spacing, .
3.1.3 Decomposing the error estimator
One application of the error estimator is to identify different contributions to the overall error, namely the error coming from discretization in space , from the discretization in time and from the splitting . This information can help to optimally balance the discretization, e.g. by avoiding excessive refinement (in space or time) or by avoiding (or by applying) a more costly IMEX (see Lemieux et. al. (2014)) approach to avoid the splitting error.
The structure of the error estimator (27) consists of residuals weighted by primal or dual interpolation errors (first and second line of (27)) and by two terms, and , which measure the splitting error . The residual terms refer to the discretization error in space and in time . An allocation of this combined spacetime error can be achieved by introducing intermediate interpolations, which we discuss for the primal residual term (the first term in (27)). Here we introduce , an interpolation into the space of functions that are piecewise constant in time but still nondiscrete in space
(29) 
We can split the residual into two separate parts, as the dependency on the weight (which takes the role of the test function) is always linear.
Naturally, the interpolation is not available. However, we can approximate the interpolation errors by the reconstruction operator that have been mentioned in Remark 3.3. To be precise, the two terms in (29) are approximated by
The philosophy is simple: for estimating the error in time, we compare the discrete solution with its higher order reconstruction in time , the spatial error is estimated by considering the spatial reconstruction operator only. All further residual terms in (27) are handled in the same way.
3.2 Realization for sea ice dynamics
In order to apply the adaptive feedback loop presented in Algorithm 3.1, we will first indicate the exact discrete formulations for solving the primal and dual problems. These are not written as dG(0) formulation but in the form of the Euler method. Afterwards we give some further remarks on the evaluation of the error estimator.
Algorithm 3.4 (Primal sea ice problem)
Let and be the initial solutions at time . Iterate for

Solve for the velocity

Solve the transport equations for and
To derive the dual sea ice model defined in Theorem 3.2, we must differentiate the form , which is described in (13), in the direction of the solution . The dual solution replaces the test function and the new test function is the argument of the directional derivative. The complete derivative of the variational formulation is presented in Mehlmann (2019).
Most characteristic for the dual problem is the reversal of the time direction, the problem runs backward in time. This property is most easily seen if we go back to the linear heat equation discussed in Section 3.1. Since the variational formulation (15) is linear, the derivative is simply the form itself and the dual form is derived by switching the role of test function and solution. We only consider the term which turns into
where we used partial integration. Up to boundary terms we obtain the dual time derivative .
This reversal of direction also carries over to the splitting scheme. While the primal iteration first solves the momentum equation, the dual problem naturally results in first solving the (dual) transport problems, followed by the (dual) momentum equation.
Algorithm 3.5 (Partitioned solution approach for the dual system)
Let and for , be the discrete solution of the primal problem. We set and and itera