Abstract

This study examines the well-posedness of the initial-value problems that arise in common models of sea ice. The model equations describe the balance of linear momentum combined with simplified thermodynamics represented by two continuity equations for effective ice thickness and ice concentration. The constitutive model for sea ice is given by two possible variants of the viscous–plastic model: the viscous–plastic model with pressure replacement and the viscous–plastic model with pressure replacement plus a tensile cutoff. The authors identify regimes of well- and ill-posedness for both models in one and two space dimensions. In one space dimension, the study finds that the viscous–plastic model and viscous–plastic model with pressure replacement behave similarly: there is ill-posedness when the divergent flow rate is larger than a minimum value. On the other hand, the viscous–plastic model with pressure replacement plus a tensile cutoff is ill-posed for all divergent flows. In two space dimensions the analysis is inconclusive for the viscous–plastic model with pressure replacement, but with the tensile cutoff the problem is ill-posed for certain divergent flows. The authors also discuss energy bounds and the difference between ill-posedness and stability of a solution. The study shows by examples that boundedness of solutions does not imply well-posedness and that it is possible for well-posed problems to have unstable solutions. The analysis shows that previous arguments in the literature, which state that a bound on the energy in sea ice models provides control over ill-posedness, are flawed.

1. Introduction

Sea ice exists at the interface between the atmosphere and ocean in the polar regions and is an important component in Earth's climate system. Current global climate models account for ice thermodynamics and ice dynamics in climate simulations. Ice dynamics is governed by a momentum balance equation. The forces acting on the ice are drag from the wind and ocean, Coriolis forces, gravitational effects from the sea surface tilt, and internal ice forces that follow a constitutive model. The constitutive model used in practically every global climate model is some variation of the viscous–plastic (VP) model introduced by Hibler (1977, 1979). This model treats ice as an isotropic, continuum, and viscous fluid at low strain rates and as a plastic solid at higher strain rates. The result of combining these behaviors is a nonlinear, viscous, and compressible fluid model for ice. In this paper, we examine an initial-value problem describing sea ice.

As with any mathematical model, it is of interest to understand if solutions to the equations have physically realistic features. Hadamard (1902) argued on physical grounds that the equations making up a mathematical model that describe a time-dependent physical process should have three features. First, a solution to the equations should exist for arbitrary initial data; second, the solution should be uniquely determined by the initial data; and third, the solution should depend continuously on the initial data. (The italicized words are vague and a precise statement of these three conditions requires specifying which class of functions is allowed and how their size is measured. We make this specification in the  appendix.) If we are attempting to describe an observed physical process with equations, then these are reasonable requirements. It is obvious that the mathematical model needs to have a solution to be useful. Uniqueness expresses the definiteness of the physical process. The third property is desirable because we usually expect that a solution based on imperfect knowledge of the initial conditions will be useful. Even if we had perfect knowledge of the initial conditions, a numerical solution of the mathematical equations necessitates approximating these conditions. Again, we would like the numerical solution to be close to the analytical solution.

If the three features introduced by Hadamard (1902) hold for an initial-value problem, then the problem is called well-posed; otherwise, it is called ill-posed. In this paper, we examine the well-posedness of an initial-value problem describing sea ice, consisting of a momentum equation plus simplified thermodynamics governed by ice concentration and thickness continuity equations. The internal forces in the momentum equation are modeled by popular variations of the VP constitutive model. Specifically, we consider the VP model with pressure replacement (Ip et al. 1991), and the VP model with pressure replacement plus a tensile cutoff (Hibler and Schulson 1997; Geiger et al. 1998). We concentrate our attention on the third property of a well-posed problem—the continuity of solutions with respect to initial data.

Well-posedness using the original VP model was studied by Gray (1999). Using a one-dimension analysis, he found that the equations with the original VP model are hyperbolic when flow is converging uniaxially, but change to mixed elliptic/hyperbolic behavior when the flow is diverging, leading to ill-posedness. We find similar behavior in one space dimension. Specifically, we find that the divergent flow rate has to be larger than a positive minimum for ill-posedness when pressure replacement is included but that all divergent flows are ill-posed when the tensile cutoff is added. Gray (1999) demonstrates a consequence of ill-posedness for numerical solutions in one space dimension by showing that solutions to the discrete equations depend strongly on the mesh size and do not converge with mesh refinement. The numerical solutions show obvious high-frequency components that cannot be controlled in the solution of the initial-value problem. This same behavior is to be expected for the models we analyze.

Section 2 presents the problems to be analyzed. The tools used in this analysis are described in section 3. We present the analysis for one space dimension in section 4 and for the full two-dimensional formulation in section 5. In two dimensions, the analysis of the linearized system of equations is based on two uncoupled systems. One system involves only the effective thickness and ice concentration. This system is well-posed. The other system involves the velocity components. With pressure replacement, but no tensile cutoff, the linear analysis of the system with velocity components is inconclusive; with the tensile cutoff, the analysis shows that the initial-value problem is ill-posed.

Questions have arisen in the past about the well-posedness of Hibler's original VP model (Hibler 1979). The difficulty in analyzing nonlinear systems of equations has led to confusion in the literature about the difference between unstable solutions and the ill-posedness of the initial-value problem. Gray and Killworth (1995) were apparently the first to notice a short wavelength linear instability in the VP model, but stated that this instability could be controlled by nonlinearities in the problem. Schulkes (1996) showed that the total energy of the system is bounded and claimed this fact made the VP model asymptotically stable. Dukowicz (1997) argued similarly that the dissipative nature of the VP model controlled the instabilities found in Gray and Killworth (1995). Gray and Killworth (1997) further responded to Dukowicz (1997) and Schulkes (1996) by correcting some flaws in the analysis and clarifying that Gray and Killworth (1995) showed ill-posedness because the growth rate of the instability increased with decreasing wavelength. But they stated erroneously that the nonlinear analysis meant that the initial-value problem was well-posed for certain boundary conditions. To make the distinction between an ill-posed initial-value problem and an unstable solution clear, we provide an  appendix defining well- and ill-posedness for the Cauchy problem with spatial differential operators of orders one and two. In section 6, we provide examples to show that (i) boundedness of solutions does not imply that the initial-value problem is well-posed and (ii) a well-posed initial-value problem can have unstable solutions. Unstable solutions can be physically meaningful; in contrast, ill-posed problems are typically not physically relevant because the models have no predictive capability.

2. Ice dynamics with viscous–plastic rheology

In this section, we present the balance equations for the prognostic variables: velocity v, ice concentration or area fraction A, and effective ice thickness h. The governing equations are

 
formula
 
formula
 
formula

where ρ > 0 is a constant ice density and σ is the depth-averaged stress. The constitutive model for σ is given below. The external forces acting on the ice are denoted by Fext. These consist of Coriolis forces, wind stress, ocean drag, and sea surface tilt. The first equation (1) describes the balance of linear momentum for a two-dimensional layer of sea ice. The reduction to two dimensions is made by integrating three-dimensional equations through the thickness. The effective thickness is the ice volume per unit area. The notation d/dt means the material time derivative d/dt = ∂/∂t + v · . The second and third equations, (2) and (3), provide a simple thermodynamic model for ice growth. These are continuity equations for the ice concentration and effective thickness, respectively. The terms SA and Sh are thermodynamic source terms. A full description of the model can be found in the references (Hibler 1979; Gray and Morland 1994). To complete the model, the differential equations are augmented by initial conditions and relevant boundary conditions.

Our goal is to investigate the well-posedness of the initial-value problem in one and two space dimensions. Therefore, we do not need a detailed description of the forcing terms Fext, SA, and Sh in (1)(3) because these terms do not contain derivatives of the dependent variables, and thus, do not affect the analysis, as we will see below.

a. The viscous–plastic constitutive law

Hibler's original model (Hibler 1979) is derivable from a rigid plastic model together with a normal flow rule. In this model, the yield curve is an ellipse in stress space given by the equation

 
formula

where σ1 and σ2 are the principal stress components of the two-dimensional stress tensor. The length of the major axis of the ellipse is in principal stress space. The ellipse is centered at (−P/2, −P/2), and e is its ratio between lengths of the major and minor axes (Fig. 1, top).

The resulting model for plastic flow (Hibler 1977) is

 
formula

in terms of the rate of deformation tensor , where DI = , and is the second-order identity tensor. (The colon notation is a contraction operation between the two second-order tensors. In components, = ijij, where repeated indices are summed.) The viscosity coefficients are

 
formula

The isotropic compressive strength P is assumed to depend on the mean ice thickness and area fraction through the relation

 
formula

with P* and C being positive constants. Usually, e = 2, C = 20, and values of P* from 2.0 × 103 to 5.0 × 104 N m−2 have been used.

The quantity Δ is a nonlinear function of the strain rate

 
formula

with DII being the maximum shear strain rate and DI being the dilation rate. Plastic flow occurs when the strain rates are large enough that the stresses remain on the elliptical curve (4) shown in Fig. 1, top. Hibler (1979) extends this formula to hold for small strain rates as well. For small strain rates, the viscosity coefficients become large because Δ gets small. Hibler (1979) cuts off Δ at a value Δmin > 0 so that the viscosity coefficients become

 
formula

where Δm = max(Δ, Δmin) with Δmin being a constant. Usually, Δmin is taken to be 2.0 × 10−9 s−1.

The model given so far, consisting of (1)(3) with the constitutive (5), is the original VP model. Well-posedness of this model in one space dimension was considered in Gray (1999) and shown to be ill-posed. However, this model has a physical flaw that is corrected by the so-called pressure replacement described in the next section.

b. Pressure replacement

At low strain rates, the viscosity coefficients do not depend on the strain rate, and the model of the previous section reduces to viscous flow. For a given strain rate in this regime, the stress lies on an ellipse in stress space that is concentric with the plastic yield curve, but smaller (Fig. 1, top). Because such an ellipse does not go through the origin, we see that we can have nonzero stress as the strain rate goes to zero. To change this nonphysical behavior, Ip et al. (1991) made the first modification to the model of the last section that is now universally accepted. The pressure is replaced with the value 2ζΔ so that the center of the ellipse is shifted (Fig. 1, middle) and the constitutive model is

 
formula

Well-posedness of this model has not been studied previously even though it is commonly implemented in large-scale climate simulations.

c. Tensile cutoff

The elliptical yield curve allows tensile stresses when −P/(e2 + 1) < σ1 < 0 or −P/(e2 + 1) < σ2 < 0. Hibler and Schulson (1997) eliminated the tensile region by reducing the shear modulus when the stress is in this region. In terms of principal stresses, we need to satisfy |σ1σ2| < |σ1 + σ2| and σ1 < 0. Therefore, to eliminate tensile stresses, it is enough to impose

 
formula

The resulting curve is illustrated in Fig. 1, bottom.

3. Overview of the analysis

We will analyze two models that differ in their definition of η.

  • Model 1 is the VP model with pressure replacement given by (1)(3) with σ given by (9), P given by (7), and ζ and η given by (8).

  • Model 2 is the VP model with pressure replacement and the tensile cutoff given by (1)(3) with σ given by (9), P given by (7), ζ given by (8), and η given by (10).

The dependent variables in these equations are the velocity v(x, y, t) with components denoted by [u(x, y, t), υ(x, y, t)], the effective thickness h(x, y, t) and the ice concentration A(x, y, t).

For either model, we introduce the vector function of unknowns

 
formula

so that both models take the general form

 
formula

where is a nonlinear, spatial differential operator, which contains derivatives up to order two for these models and F contains the forcing terms. We are interested in the pure initial-value problem, or Cauchy problem, consisting of the system of differential equations with prescribed initial conditions

 
formula
 
formula

for x, y ∈ ℝ, and for some interval in time 0 ≤ tT. We also assume that this problem has a unique solution denoted by U0(x, y, t) = [u0(x, y, t), υ0(x, y, t), A0(x, y, t), h0(x, y, t)]T for h0 > 0.

For the problem in (11) and (12) to be well-posed at U0, we ask that this solution depend continuously on the data F and f. That is, if we perturb the data by replacing F and f with and , then the perturbed problem is also uniquely solvable, at least for small perturbations. Moreover, if we write the perturbed solution as

 
formula

then the perturbation in the solution should be bounded by a constant times the perturbation in the data, in appropriate norms. Note that different norms lead to different notions of well-posedness, as discussed in the  appendix. The question of well-posedness of the nonlinear problem is closely related to well-posedness of the linearized equations. If all linear problems obtained by linearizing the nonlinear problem at all functions near U0 are well-posed then, under certain circumstances, the nonlinear problem is also well-posed at U0. It also suffices to consider the homogeneous problem given by F = 0 because the forcing terms do not contain derivatives of the dependent variables. Thus, these lower-order terms do not contribute to the analysis.

The linearized equations for the homogeneous problem are obtained by substituting (13) into (11), with F = 0, and formally neglecting terms of order ɛ2, to obtain a linear set of equations for . For our models, the general form of the linearized system of equations is

 
formula

where the coefficient matrices ij, j, and depend on the solution U0 and its derivatives.

We note that the operator is a linear differential operator and (14) is a linear equation with nonconstant coefficients. The analysis proceeds by examining this problem when we freeze the coefficients, a process called localization. The motivation is the following: if all problems with frozen coefficients are well-posed then, under some additional assumptions, the corresponding variable coefficient problem is also well-posed. This localization principle is not true in general, but for certain classes of operators it is true (Kreiss and Lorenz 2004). The  appendix summarizes results on well-posedness for parabolic, hyperbolic, and mixed hyperbolic–parabolic-type problems that we will use to analyze our model problems. For these operators, the localization principle is valid (Kreiss and Lorenz 2004).

4. One-dimensional flow

We begin the analysis by considering one-dimensional problems. Of course, the same process described in the last section can be applied to uniaxial flow in one space dimension. In that case, the vector function U does not depend on y and the variable υ does not appear. In one dimension, the homogeneous form of (1)(3) is

 
formula
 
formula
 
formula

The quantity σ is the xx component of the stress tensor, σ = (ζ + η)uxζΔ. In one space dimension, we have the following simplifications: D = ux, DI = ux, DII = |ux|, and , where E = 1 + 1/e2.

a. Model 1

The first model is the VP model with pressure replacement. In this case, substitute ζ and η from (8) to get the xx component of the stress

 
formula

Because Δm = max(Δmin, Δ), we next consider the two cases (i) Δ < Δmin and (ii) Δ > Δmin, separately. In each case, we assume that the base flow U0(x, t) falls in the indicated regime and that the perturbation is small enough that the perturbed flow, , remains within the same regime.

1) Case Δ < Δmin

In one dimension, the inequality defining the case under consideration reduces to and the xx component of stress can be written as σ = SPux, where

 
formula

Following the prescription of the last section, we want to study the linear system of equations for the perturbation . We assume that the perturbation is small enough that the sign of ux is the same as the sign of u0x. Dropping the tildes, the linearized system of equations becomes

 
formula

where the linearized pressure perturbation is P = P0(h/h0 + ) with . This system is a coupled hyperbolic–parabolic system, which is well-posed according to the  appendix (section d). The uncoupled hyperbolic part is

 
formula

and the parabolic part is ut = (SP0/ρh0)uxx, with SP0/(ρh0) > 0.

In one space dimension, the regime under consideration, Δ < Δmin, reduces to . Thus, in the viscous flow regime where , the VP model with pressure replacement leads to a well-posed system of equations.

2) Case Δ > Δmin

In this case, the xx component of the stress is . We assume that the perturbation is small enough so that the sign of ux is the same as the sign of u0x. The linearized system of equations becomes

 
formula

where, as before, P = P0(h/h0 + A) with , and we have dropped the tildes from the perturbed quantities to simplify the notation.

The linear system has the form Ut = 1Ux + U where U = (u, A, h)T. In this case, set

 
formula

and write out the matrix

 
formula

The characteristic equation of this matrix is

 
formula

Note that S ≠ 0 and P0(1 + A0)/ρ > 0. When S < 0, the characteristic equation has three distinct real roots. However, when S > 0, the characteristic equation has two roots with a nontrivial imaginary part. According to the  appendix (section b), the Cauchy problem is well-posed when S < 0 because the equations are strictly hyperbolic. However, the problem is ill-posed when S > 0. Notice that the flow is converging if S < 0 and the flow is diverging when S > 0. Our result is the same as the result in Gray (1999) for the original VP model; with pressure replacement, the problem is ill-posed when the divergent flow rate is larger than a minimum value, .

The results for model 1 are summarized in Table 1. The physical mechanisms that lead to ill-posedness in divergent flow, discussed in Gray (1999) and Gray and Killworth (1995), remain valid. There is a positive feedback loop created by coupling between the time derivative of the velocity and the strain hardening constitutive model through the ice concentration and thickness equations. If there is a weak area, then P is less, there is less resistance to motion, the ice diverges more, and A and h decrease further, establishing the feedback makes the small weakness grow.

The advection terms are not important in this feedback loop. If we examine the dominant terms given by Ut1Ux = 0, but drop the advection terms (or examine the equations in a frame moving with the flow), then we have the following equations

 
formula

We can eliminate A and h by differentiating the first equation with respect to t and the last two equations with respect to x. The result is a single equation for the velocity u

 
formula

When S < 0, this is a wave equation for the velocity and an initial-value problem makes sense. When S > 0, this is an elliptic problem for the velocity, like Laplace's equation. It does not make sense to solve Laplace's equation as an initial-value problem and, physically, we do not expect the velocity to be governed by such an equation. This observation is the essence of Hadamard's argument that ill-posed problems are malformed.

b. Model 2

In this section, we examine the VP model with pressure replacement plus the tensile cutoff. In one dimension, η given in (10) simplifies to

 
formula

with ζ = P/(2Δm) and Δm = min(Δ, Δmin). It is convenient to divide the flow into four regimes: regime 1 , regime 2 , regime 3 , and regime 4 . In regimes 1 and 2, ux < 0 and η and ζ are the same as in model 1. The analysis of case 1 in the last section shows that the problem of this section is well-posed in regime 2; and the analysis of case 2 in the last section shows that the problem of this section is well-posed in regime 1. The tensile cutoff changes the analysis for regimes 3 and 4. In these regimes, ux > 0, , and therefore, σ = (ζ + η)uxζΔ = 0, as designed by the tensile cutoff. Although the definition of ζ is different in regimes 3 and 4, in one dimension the value of ζ plays no role in these regimes and they can be analyzed together.

We assume that the perturbation is small enough that ux falls in the same regime as u0x. The linearized system of equations becomes

 
formula

where, as before, we have dropped the tildes from the perturbed quantities to simplify the notation. This system has the form Ut = 1Ux + U where the matrix 1 is

 
formula

The Jordan form (e.g., Strang 2007) of −1 is

 
formula

As stated in the  appendix (section b), the presence of a 2 × 2 Jordan block implies ill-posedness according to definition A1 for the constant coefficient problem. Regimes 3 and 4 correspond to diverging flow. The ill-posedness for divergent flow is more subtle for model 2 than model 1. For model 2, the dominant terms given by Ut = 1Ux constitute a system that is weakly well-posed in the sense of (A10). The weak well-posedness and the presence of lower-order terms in (20) then makes the nonlinear problem ill-posed (in any sense). The growth of perturbations is caused by the term involving h in the velocity equation.

The results of this section are summarized in Table 1. The tensile cutoff does not remedy the ill-posedness in the VP model with pressure replacement. In fact, the condition for ill-posedness becomes ux > 0.

5. Two-dimensional flow

In this section, we analyze the full two-dimensional equations given in section 2 with the constitutive models denoted above as models 1 and 2. The first is the VP model with pressure replacement and the second includes the tensile cutoff. In either case, the linearized equations can be written

 
formula

where

 
formula

The terms in the last equation for the perturbed stress have the following definitions:

 
formula

We also note the pressure expansion

 
formula

with . The exact form of and depends on the specific model under consideration and will be given below. As stated previously, the exact form of the external forces, , , , is not required.

There are no second derivatives of and in these equations. All second-order spatial derivatives arise from differentiation of and appear as second derivatives of . Therefore, the system of linearized equations has the following structure:

 
formula

where we have dropped the tildes. The matrices have a block structure with each block being 2 × 2. The matrix is the 2 × 2 zero matrix.

In the remainder of this section, we examine whether or not this system of equations constitutes a coupled hyperbolic–parabolic system as defined in the  appendix (section g). To this end, we first consider

 
formula

According to definition A4, it is easy to check if this system is strongly hyperbolic. Next we must check if

 
formula

is parabolic according to definition A5. This condition is checked for each model separately in the following.

a. Model 1

In this section we consider the VP model with pressure replacement. As in one dimension, the two cases (i) Δ < Δmin and (ii) Δ > Δmin are analyzed separately.

1) Case Δ < Δmin

In this case, ζ and η take the form

 
formula
 
formula

The matrices that appear in (24) are

 
formula

with E = 1 + 1/e2, as before.

According to definition A5, parabolicity of (24) requires that the eigenvalues λl of the matrix satisfy

It suffices to test this condition for |ω| = 1, ω = (ω1, ω2), because the eigenproblem can be scaled by |ω|. We test this condition numerically.

We let each of u0x, υ0x, u0y, and υ0y take values from the set (−3s, −2s, −1s, 0, 1s, 2s, 3s), with s serving as a scaling factor and taking values s = 10−8, 10−9, 10−10, and 10−11. In other words, we sample 74 = 2401 values of the vector (u0x, υ0x, u0y, υ0y). For each of the vectors (u0x, υ0x, u0y, υ0y), we verify Δ0 < Δmin = 2.0 × 10−9. If the inequality holds, we mark the vector (u0x, υ0x, u0y, υ0y) as qualified and obtain the matrices jk. Next, we form a matrix for each of 50 pairs of ω1 = cos(θ), ω2 = sin(θ), where θ varies from 0 to 2π in steps of 2π/50. The real parts of the eigenvalues of are recorded. Finally, the minimum of the real part of the eigenvalues for all possible qualified values of (u0x, υ0x, u0y, υ0y) and each (ω1, ω2) is obtained. This computation leads to the conclusion

 
formula

The numerical evidence suggests that the system (22) with frozen coefficients is well-posed in this case; and therefore the VP model with pressure replacement and Δ < Δmin is well-posed.

2) Case Δ > Δmin

In this case, ζ and η take the form

 
formula
 
formula

The matrices that appear in (24) are

 
formula

With the one-dimensional analysis in mind, consider a uniaxial flow with u0x ≠ 0 and u0y = υ0x = υ0y = 0. Then we obtain DI = u0x, DII = |u0x|, , and

 
formula

The matrix

 
formula

has a zero eigenvalue when ω1 = 1 and ω2 = 0.

Also, we explore the parabolicity condition numerically, assuming that each of u0x, υ0x, u0y, and υ0y takes values from (−3s, −2s, −1s, 0, 1s, 2s, 3s), and varying ω as before. Considering s = 10−3, 10−4, 10−8, and 10−9 and taking into account only values of u0x, υ0x, u0y, and υ0y that obey Δ0 > Δmin = 2.0 × 10−9, we obtain that eigenvalues λl of the matrix satisfy up to round-off error. We have not been able to find parameter values that yield eigenvalues with a negative real part. Thus, our analysis is inconclusive in this case.

b. Model 2

In this section, we consider the VP model with pressure replacement and the tensile cutoff. Recalling that , we see that (Δ − DI)/DII > 1/e2 if and only if DI < 0. Thus, we are led to consider four flow regimes:

  • Regime 1: Δ > Δmin and DI < 0, (η = ζ/e2)

  • Regime 2: Δ < Δmin and DI < 0, (η = ζ/e2)

  • Regime 3: Δ < Δmin and DI ≥ 0, [η = ζ(Δ − DI)/DII]

  • Regime 4: Δ > Δmin and DI ≥ 0, [η = ζ(Δ − DI)/DII]

Regimes 1 and 2 use the same values of ζ and η as model 1. Uniaxial flow with u0x < 0 falls into regime 1 and we can expect a zero eigenvalue of the matrix . Regime 2 is well-posed. For regimes 3 and 4, we again explore the test for parabolicity numerically. Because of tedious algebra, for these regimes, we calculate the derivatives in the matrices jk using central finite differences of second order. We consider each of u0x, υ0x, u0y, and υ0y, taking values from (−3s, −2s, −1s, 0, 1s, 2s, 3s), and ω1 = cos(θ), ω2 = sin(θ), where θ ∈ [0, 2π) in steps of 2π/50. For regime 3, we use s = 10−8, 10−9, 10−10, and 10−11. The qualified values of u0x, υ0x, u0y, and υ0y are those obeying Δ0 < Δmin = 2.0 × 10−9. We obtain min(Reλ1,2) ≃ 0.5. Similarly, for regime 4 we use the same values for the variables but with s = 10−2, 10−3, and 10−9. Qualified values of u0x, υ0x, u0y, and υ0y obey the criterion Δ0 > Δmin = 2.0 × 10−9. We obtain min(Reλ1,2) ≃ −23.5. Thus, the numerical evidence suggests the model is well-posed in regime 3, but ill-posed for flows in regime 4. The results of this section are summarized in Table 2.

6. Comments on dissipation, boundedness of solutions, and stability

Ill-posed problems are not physically realistic (Hadamard 1902). On the other hand, there are many examples of physically relevant unstable flows (e.g., Drazin and Reid 1981). The fact that perturbed solutions are unstable and therefore grow in time is not equivalent to a problem being ill-posed. Establishing well- or ill-posedness requires more than just examining growth or decay of perturbations—it requires the presence or absence of the bound (A7) in definition A1. The analysis of stability of solutions begins similarly to the analysis of well-posedness. Given a solution U0, perturb it by , obtain a linearized equation for the perturbation, and ask whether or not the perturbation grows or decays in time. To illustrate the process, consider the following model problem in one dimension,

 
formula

Now consider a solution u0(x, t) to the initial-value problem (29) and perturb the initial data The solution to the perturbed problem is denoted . The perturbation solves the initial-value problem with . If the initial perturbation is a single Fourier mode, , then . The perturbation grows unboundedly in time if |ω| < 1. More general perturbations can be handled using Fourier transformation and the conclusion is the same, u0 is unstable if the perturbation contains modes with |ω| < 1.

However, this problem is well-posed at u0. For this example, the symbol of the operator is . Because , this problem is well-posed according to definition A1. We see that it is possible for solutions to a well-posed problem to have perturbations that grow unboundedly in time; that is, for the solution u0 to be unstable. Because of the bound (A7), well-posedness requires that growth not occur for modes with |ω| → ∞; such growth allows solutions to blow up in finite time for the linearized equation. Moreover, the presence or the absence of the bound (A7) allows, under circumstances noted in the  appendix, for conclusions about well-posedness of the linearized, frozen coefficient problem to be extended to the original nonlinear problem. It is this property that allows us to draw conclusions about the sea ice models analyzed in sections 4 and 5.

A dissipative constitutive model allows one to show that the energy in the system is bounded. This property is also desirable for modeling physically realistic systems. A constitutive model that is not dissipative violates the second law of thermodynamics and should be discarded on physical grounds. Dukowicz (1997) showed that the original VP model is dissipative during plastic flow. However, Gray and Killworth (1997) and Schulkes (1996) showed the model is not dissipative in the viscous regime for divergent flow. Pritchard (2005) showed that the VP model with pressure replacement is always dissipative but the VP model with pressure replacement plus the tensile cutoff is not dissipative if DI > DII. Hibler and Schulson (2000) also noted this fact.

Dukowicz (1997), Schulkes (1996), and Gray and Killworth (1997) argue that a bound on the energy provides control over the ill-posedness of the sea ice model. Their argument is that the linear analysis governs the initial growth of perturbations but once these perturbations grow, nonlinearities in the original problem must control that growth because the energy is bounded. This argument is flawed. We next give an example of a nonlinear, ill-posed initial-value problem where an energy bound can be established. The example demonstrates that an energy bound does not suffice to show well-posedness. Consider the one-dimensional initial-boundary-value problem

 
formula

where f(x) is a smooth, real-value function satisfying the boundary conditions. We will first show that the energy of every smooth solution remains bounded. Nevertheless, by considering initial data

 
formula

where c is a small positive constant and j is large, we will show that there is no bound on the growth rate of the solutions, and the problem is ill-posed.

Using integration by parts, we have for the change of energy

 
formula

For any fixed t and at any point x, the integrand of the integral J does not exceed the value one. Indeed, the integrand is a polynomial in ux, namely, . This polynomial is zero when ux = 0 or 1, and has a maximum value of 3/16 when ux = 1/2. Therefore, J ≤ 1, and we have shown that d(t)/dt ≤ 1 − (t). For the function h(t) = et(t) we obtain that h′(t) = et′(t) + et(t) ≤ et, thus h(t) ≤ (0) + et. This proves the bound (t) = eth(t) ≤ et(0) + 1 for all t ≥ 0.

To show ill-posedness of the problem, we write the equation in the form and consider the initial condition (30) where c = 1/3π and j is large. We have ux(x, 0) = cos(jπx), thus It follows that the initial-value problem is as badly behaved as the backward heat equation, at least for a short time and small initial data. As j → ∞, there is no bound for the growth rate of solutions with initial data (30).

The one-dimensional numerical examples in Gray (1999) show the combined effects of ill-posedness and an energy bound for ice dynamics with the original VP rheology. Round-off errors introduce perturbations of the concentration. The perturbations grow in amplitude with increasing time, appearing as narrow fingers of decreasing ice concentration. The spacing between the fingers depends on the mesh size because the highest-frequency perturbations grow fastest according to the analysis. Thus, numerical solutions depend on the mesh size and the notion of convergence of solutions with mesh size is lost. The energy bound means that these perturbations do not grow forever and solutions remain finite. The same features are expected for the VP model with pressure replacement in one space dimension. Because the VP model with pressure replacement plus the tensile cutoff is not dissipative, the long-term growth of perturbations can potentially be more severe.

7. Conclusions

We have examined well-posedness of two models. Model 1 is the VP model with pressure replacement and model 2 is the VP model with pressure replacement plus the tensile cutoff. In one space dimension, we are able to give a complete analysis. Model 1 is ill-posed when , that is, for divergent flow when the flow rate is above a minimum value. This finding agrees with analysis in Gray (1999) for the original VP model. Model 2 is ill-posed when ux > 0, that is, for all divergent flow. The results for one space dimension are summarized in Table 1. In two space dimensions, the analysis consists of examining an uncoupled system of equations for A and h, and an uncoupled system for the velocity components. The uncoupled equations for A and h are always strongly hyperbolic and therefore well-posed. Thus, well-posedness hinges on the analysis of the uncoupled system for the velocity components. The analysis of the velocity equations is more difficult and we rely mainly on numerical testing of the requisite conditions. In two space dimensions, the numerical evidence suggests that model 1 is well-posed when Δ < Δmin. The numerical analysis of model 1 is inconclusive when Δ > Δmin. However, if in regime 4, where Δ > Δmin and DI > 0, the flow is one-dimensional, then the problem is ill-posed based on the analysis in one dimension. On the other hand, for model 2, the numerical evidence suggests the problem is ill-posed whenever Δ > Δmin and DI > 0, that is, for divergent flow. Table 2 summarizes the results for two space dimensions.

In this paper, we have tried to make clear the distinction between unstable flows and an ill-posed problem. In particular, we showed by an example that it is possible for an ill-posed problem to have a bound on the solution energy. Moreover, it is possible to have a well-posed problem with unstable solutions. Unstable solutions can be physically meaningful but the explosive, uncontrollable growth associated with perturbations of an ill-posed problem is generally not physical. Ill-posedness means that the behavior of solutions to the initial-value problem can be arbitrarily bad, and the model has no predictive capability. Because numerical approximations introduce round-off errors, which can be thought of as high-frequency perturbations to the solution, standard numerical methods cannot be expected to converge for such models. These properties are clearly undesirable. A particularly insidious issue for an ill-posed sea ice model with bounded energy is that numerical solutions can seem reasonable. The exponential growth of perturbations to the linearized system does not appear in the numerical solution of the nonlinear problem because it is quenched by the nonlinear terms responsible for the energy bound. However, the ill-posedness still indicates a badly formed problem and numerical solutions are meaningless at any resolution. For the VP model, the mechanism leading to uncontrolled growth and the implications of bounded energy are identified in Gray (1999).

From a mathematical standpoint, it is possible to “regularize” ill-posed models to make them well-posed. One approach is to add dissipative terms to the equations to kill off high-frequency growing modes. However, one has to be careful that the regularized model still represents correct physics. Hibler (1979) added harmonic and biharmonic diffusion terms to the thickness and ice concentration equations, which he stated “is viewed primarily as a necessary numerical artifact to suppress smaller scales” (Hibler 1979, p. 822). The diffusion coefficients are taken to be dependent on the grid size. In this case, it is not clear what the limit of the discrete equations is as the mesh is refined. Because well-posedness is a property of the differential equations, the differential equations need to be regularized and then appropriate numerical methods can be designed to converge to the regularized equations. Gray (1999) suggests changing the model to include rate effects as a regularization, or by shifting the yield curve. These modifications are still viable.

Acknowledgments

This material is based upon work supported by the National Science Foundation under Grant ARC-1023667.

APPENDIX

Well-Posedness

a. General concepts

In this appendix, we summarize definitions and results on well-posedness used in the paper. We refer to Kreiss and Lorenz (2004) for details and proofs.

We specialize the general theory to spatial differential operators of order one or two and consider the Cauchy problem for a system of PDEs:

 
formula

where is a constant coefficient operator of the general form

 
formula

Here, U = U(x, t) is a vector of length n; jk, j, ∈ ℂn×n are constant matrices; and Dj = ∂/∂xj.

Fix ω ∈ ℝs and consider an initial function of the form

 
formula

where is a fixed vector and .To obtain a solution of (A1) with initial data (A3), use the ansatz

 
formula

and find that

 
formula

where

 
formula

is the so-called symbol of the operator . The symbol is obtained from by formally substituting j for Dj. It then follows that

 
formula

solves the Cauchy problem (A1) with initial data (A3). Solutions of (A1) with more general initial data f(x) can be obtained by Fourier transformation. Together with Parseval's relation, this motivates the following definition.

Definition A1: The Cauchy problem (A1) is called well-posed if there exist real constants K and α so that

 
formula

This definition uses the standard notion of a vector norm |·| and its corresponding matrix norm. Specifically, for U ∈ ℂn define the norm as |U| = 〈U, U1/2, where the inner product of any two complex-valued vectors, U, V ∈ ℂn, is given in the standard way by . The corresponding norm of a matrix ∈ ℂn×n is defined as || = max(|U|, |U| = 1).

Remarks:

  1. Well-posedness according to definition A1 is equivalent to the solution estimate 
    formula
    where 
    formula
    denotes the spatial L2 norm. Other concepts of well-posedness are also in use. For example, replacing (A7) by the weaker bound 
    formula
    amounts to solution estimates of the form 
    formula
    where first derivatives of the initial function f(x) = U(x, 0) are needed to bound the spatial L2 norm of the solution U(x, t) at a later time t > 0. It turns out that the well-posedness concept of definition A1 has advantages, however: It is possible to deduce well-posedness of variable coefficient problems—and even suitable nonlinear problems—from the well-posedness of corresponding frozen coefficient problems. This is not possible, in general, if one allows weaker bounds like (A10) for frozen coefficient problems. We refer to Kreiss and Lorenz (2004) for details.
  2. The theory of well-posedness for PDEs becomes more difficult if one considers initial-boundary-value problems instead of the Cauchy problem. Even if the Cauchy problem for a system Ut = U is well-posed, the introduction of boundary conditions may lead to ill-posedness. However, if the Cauchy problem is ill-posed, the introduction of boundary conditions typically does not lead to a well-posed initial-boundary-value problem. This is plausible for the following reason: ill-posedness of the Cauchy problem means that the matrix exponential does not satisfy a bound (A7) for large |ω|. Because large values of |ω| correspond to high-frequency functions (A3), the explosion in time is of a local nature in space. The explosion cannot be prevented by (far away) boundaries.

b. One space dimension, first order

Consider a system

 
formula

where 1 and are constant matrices. The following terminology is standard:

Definition A2: The first-order system Ut = 1Ux + U is called

  • weakly hyperbolic if the eigenvalues of 1 are real;

  • strictly hyperbolic if the eigenvalues of 1 are real and distinct; and

  • strongly hyperbolic if the eigenvalues of 1 are real and 1 has a complete set of eigenvectors, that is, 1 is similar to a real diagonal matrix.

Clearly, strict implies strong hyperbolicity. One can prove that strong hyperbolicity of the system Ut = 1Ux + U is equivalent to well-posedness of the Cauchy problem (A11) according to definition A1. If 1 has a nonreal eigenvalue then (A11) is ill-posed. Also, if all eigenvalues of 1 are real but 1 is not diagonalizable, then (A11) is ill-posed. Weak well-posedness with a 2 × 2 Jordan block of the matrix 1 leads to an estimate (A10) which, generally, implies ill-posedness according to definition A1.

c. One space dimension, second order

Consider a system

 
formula

where 2, 1, and are constant matrices.

Definition A3: The second-order system Ut = 2Uxx + 1Ux + U is called parabolic if the eigenvalues λj of 2 satisfy

 
formula

One can show that the Cauchy problem for a parabolic system is well-posed according to definition A1. If 2 has an eigenvalue λj with Reλj < 0, then the Cauchy problem (A12) is ill-posed. If the eigenvalues of 2 satisfy Reλj ≥ 0, then (A13) may be well- or ill-posed.

d. One space dimension, coupled hyperbolic–parabolic systems

Consider a differential system for a vector function

 
formula

of the form

 
formula

with constant coefficient matrices. If the system Ut = 2Uxx is parabolic and the system Vt = 1Vx is strongly hyperbolic, then the Cauchy problem for the coupled system (A14) is well-posed.

e. Two space dimensions, first order

Consider a system

 
formula

where 1, 2, are constant matrices. Well-posedness does not depend on the matrix and we will assume = 0. Then the symbol of is the matrix

Definition A4: The first-order system Ut = 1Ux + 2Uy is called strongly hyperbolic if the following two conditions hold:

  1. For all (ω1, ω2) ∈ ℝ2, the eigenvalues of ω11 + ω22 are real.

  2. There is a constant K and for all (ω1, ω2) ∈ ℝ2 there is a transformation S(ω) ∈ ℂn×n so that is diagonal and |S(ω)| + |S−1(ω)| ≤ K.

One can show that the Cauchy problem (A15) is well-posed according to definition A1 if and only if the system Ut = U is strongly hyperbolic.

A simple situation occurs if (A15) is scalar. Then, the second condition in definition A4 is trivially satisfied. Consequently, if (A15) is scalar, then the equation is strongly hyperbolic if and only if 1 and 2 are real.

f. Two space dimensions, second order

Consider a system

 
formula

where the matrices jk, j, and are constant.

Definition A5: The above second-order system Ut = U is called parabolic if there exists δ > 0 so that the eigenvalues λl of the matrix satisfy Reλlδ|ω|2 for all ω ∈ ℝ2.

Thus, parabolicity of the system Ut = U is independent of the lower-order term 1Ux + 2Uy + U.

One can show that the Cauchy problem for a parabolic system (A16) is well-posed according to definition A1.

If, for some ω ∈ ℝ2 the matrix has an eigenvalue λl with negative real part, then the Cauchy problem for Ut = U is ill-posed. In the borderline case, where eigenvalues with zero real parts occur, the Cauchy problem may be ill- or well-posed.

g. Two space dimensions, coupled hyperbolic–parabolic systems

Consider a coupled system, with constant coefficients, of the form

 
formula

where Ut = 2U is second-order parabolic; Vt = 1V is first-order strongly hyperbolic; and R1 and R2 are first-order operators. The Cauchy problem for such a coupled hyperbolic–parabolic system is well-posed.

REFERENCES

REFERENCES
Drazin
,
P.
, and
W.
Reid
,
1981
:
Hydrodynamic Stability.
Cambridge University Press, 628 pp.
Dukowicz
,
J.
,
1997
:
Comments on “Stability of the viscous–plastic sea ice rheology.”
J. Phys. Oceanogr.
,
27
,
480
481
.
Geiger
,
C. A.
,
W. D.
Hibler
III
, and
S. F.
Ackley
,
1998
:
Large-scale sea ice drift and deformation: Comparison between models and observations in the western Weddell Sea during 1993
.
J. Geophys. Res.
,
103
(
C10
),
21 893
21 913
.
Gray
,
J.
,
1999
:
Loss of hyperbolicity and ill-posedness of the viscous–plastic sea ice rheology in uniaxial divergent flow
.
J. Phys. Oceanogr.
,
29
,
2920
2929
.
Gray
,
J.
, and
L.
Morland
,
1994
:
A two-dimensional model for the dynamics of sea-ice
.
Phys. Eng. Sci.
,
347
,
219
290
.
Gray
,
J.
, and
P.
Killworth
,
1995
:
Stability of the viscous–plastic sea ice rheology
.
J. Phys. Oceanogr.
,
25
,
971
978
.
Gray
,
J.
, and
P.
Killworth
,
1997
:
Reply
.
J. Phys. Oceanogr.
,
27
,
482
483
.
Hadamard
,
J.
,
1902
:
Sur les problèmes aux derivées partielles et leur signification physique
.
Princeton Univ. Bull.
,
XIII
,
49
52
.
Hibler
,
W.
,
1977
:
Viscous sea ice law as a stochastic average of plasticity
.
J. Geophys. Res.
,
82
,
3932
3938
.
Hibler
,
W.
,
1979
:
Dynamic thermodynamic sea ice model
.
J. Phys. Oceanogr.
,
9
,
815
846
.
Hibler
,
W.
, and
E.
Schulson
,
1997
:
On modeling sea-ice fracture and flow in numerical investigations of climate
.
Ann. Glaciol.
,
25
,
26
32
.
Hibler
,
W.
, and
E.
Schulson
,
2000
:
On modeling the anisotropic failure and flow of flawed sea ice
.
J. Geophys. Res.
,
105
(
C7
),
17 105
17 120
.
Ip
,
C.
,
W.
Hibler
, and
G.
Flato
,
1991
:
On the effect of rheology on seasonal sea-ice simulations
.
Ann. Glaciol.
,
15
,
17
25
.
Kreiss
,
H.-O.
, and
J.
Lorenz
,
2004
:
Initial-boundary value problems and the Navier–Stokes equations. Classics in Applied Mathematics, Vol. 47. Society for Industrial and Applied Mathematics (SIAM), 402 pp. + xviii
.
Pritchard
,
R.
,
2005
: Stability of sea ice dynamics models: Viscous–plastic rheology, replacement closure, and tensile cutoff.
J. Geophys. Res.
, 110, C12010,
doi:10.1029/2003JC001875
.
Schulkes
,
R.
,
1996
:
Asymptotic stability of the viscous–plastic sea ice rheology
.
J. Phys. Oceanogr.
,
26
,
279
283
.
Strang
,
G.
,
2007
: Linear Algebra and Its Applications: International Edition. 4th ed. Cengage Learning, 576 pp.

Footnotes

*

Sandia National Laboratories is a multiprogram laboratory operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy's National Nuclear Security Administration under Contract DE-AC04-94AL85000.