Landfast ice is sea ice that forms and remains fixed along a coast, where it is either attached to the shore or held between shoals or grounded icebergs. The current generation of sea ice models is not capable of reproducing certain aspects of landfast ice behavior, for example the persistence of landfast sea ice under the effect of offshore winds. The authors develop a landfast sea ice model by adding tensile strength to the viscous–plastic as well as two versions of the elastic–viscous–plastic sea ice rheologies. One-dimensional implementations of these rheologies are used to explore the ability of coastal sea ice to resist offshore winds over extended times. While all modified rheologies are capable of maintaining landfast ice–like structures in the model, only the viscous–plastic rheology fulfills theoretical expectations. The elastic–viscous–plastic rheologies show initial elastic waves that weaken the ice and thus reduce its capacity of maintaining landfast ice. Further, special care has to be taken when implementing the most commonly used version of the elastic–viscous–plastic rheology because the standard set of parameters is not adequate for landfast sea ice modeling.
Landfast ice is sea ice that forms and remains fixed along a coast, where it is either attached to the shore or held between shoals or grounded icebergs. It covers Arctic shelves during large portions of the year, normally starting to form in October and reaching its widest extent during April and May (Volkov et al. 2002; Barber and Hanesiak 2004; Divine et al. 2005). Because of its lack of mobility, landfast ice fundamentally modifies the momentum exchange between atmosphere and ocean. It covers the shelf area for four to eight months each year and thus greatly affects the heat and freshwater exchange there. The freezing and melting of landfast ice make important contributions to salt and freshwater budgets, thereby influencing water circulation, dense water production, and the location of upwelling and downwelling zones (Macdonald et al. 1999). In addition to its strong influence on the coastal oceanic ecosystem, the landfast ice edge is also essential to the subsistence of numerous Arctic mammals and Inupiat (Tynan and DeMaster 1997).
While available landfast sea ice data are very sparse, several claims have been made over time with respect to where it can be found. Wadhams (1986) states that the Arctic landfast sea ice edge will follow the 18-m isobath with occasional bulges out to 27 m, whereas according to Barber and Hanesiak (2004) during the northern winter it will approach the shelf break region at about 20-m ocean depth. Divine et al. (2004, 2005) find in the Kara Sea that fast ice extends to the 10-, 20–30-, or even 100-m bathymetry line, depending on the region. Mahoney et al. (2007a) state that the Arctic landfast sea ice edge will roughly follow the 18-m isobath in midseason and reach the 16–19-m water depth toward the end of the season, based on their landfast sea ice climatology for the years 1996–2004; however, they also point out that the landfast sea ice edge can reach into water depths of 3500 m. König Beatty (2007) provides a landfast sea ice climatology using data from 1953–98 in which it can be seen that fast ice in the Arctic often appears on depths much deeper than 20 m, particularly between islands of the Canadian Archipelago and on Russian continental shelves. Fast ice observations and data in the Southern Hemisphere are even sparser, the fast-ice maps for 1997 and 1999 by Giles et al. (2008) being the first for East Antarctica.
Among all the conflicting claims of where landfast sea ice can be found, there is never a reason given why the fast ice edge should follow a certain depth, and isobaths are exclusively used as a simplified description. While observations indicate that landfast sea ice seems to prefer shallow water, there are many exceptions to that rule. Hence, bathymetry alone is insufficient for accurate prediction of landfast sea ice occurrence.
Two mechanisms are generally thought to accelerate the formation of landfast sea ice, while three mechanisms seem to be responsible for its maintenance. First, initial landfast sea ice formation is accelerated by shallow ocean depths that allow faster local freezing of water due to the lack of deeper warmer water acting as a source of heat through convection. The second accelerating mechanism is integration and consolidation of pack ice that is being imported by onshore winds (Wadhams 1986). Once established, grounded ice ridges stabilize the landfast sea ice. Additionally, restrictive geometry (like groups of islands, narrow passages, and concave coastlines) enables sea ice to remain fast. The first mechanism is particularly important in shallow areas up to about 20-m to at most 25-m ocean depths (e.g., Mahoney et al. 2007b, for Point Barrow, Alaska). The second seems dominant in the Kara Sea where groups of islands stabilize the fast ice (Divine et al. 2004, 2005), in the narrow straits of the Canadian Arctic Archipelago, along Greenland, and in much of Antarctica (Wadhams 1986; Heil 2006). Tremblay and Hakakian (2006) actually find that sea ice exhibits tensile strength when subjected to offshore winds, particularly in the East Siberian Sea. This suggests a third mechanism for landfast ice maintenance.
To improve large-scale modeling of landfast sea ice, these mechanisms or an appropriate parameterization thereof have to be included in the sea ice description. The coarse spatial resolution of basinwide or global models will not allow a realistic description of the grounding of ice keels any time soon. Although bathymetry and coastline geometry are already included in sea ice models, current sea ice models fail to create realistic landfast sea ice behavior; the shear strengths used are not sufficient for fast ice to remain attached to a coastline under the influence of offshore winds, even though this phenomenon is observed in the field (Thorndike and Colony 1982; Divine et al. 2005; Heil 2006; Tremblay and Hakakian 2006).
Previous modeling attempts have concentrated on one-dimensional (in the vertical) modeling of the thermodynamics of landfast sea ice in the Arctic (Flato and Brown 1996; Dumas et al. 2006) and the Antarctic (Heil et al. 1996). Lieser (2004) implemented a scheme in which sea ice became “fast” and thus immobile whenever the sea ice thickness would reach or exceed 10% of the local ocean depth up to a maximum ocean depth of 30 m. To our knowledge no dynamical landfast sea ice modeling has been attempted so far.
Commonly, it is assumed that sea ice has no or a very low tensile strength (e.g., Hibler 1979) as the ice is finely fractured on scales that are typically used for sea ice models. Although Zhang and Rothrock (2005) studied the effect of tensile strength, the amounts they used were small and not suitable for modeling landfast sea ice. Adding tensile strength can be thought of as a parameterization for subgrid-scale grounding of ice ridges, at least in areas of shallow water. Recent analysis of data by Tremblay and Hakakian (2006) and by König Beatty (2007) has shown that landfast sea ice can resist strong offshore winds, even in regions where the ocean is too deep for shoals or grounded icebergs to play a stabilizing role. Thus, the assumption of no tensile strength seems flawed in the case of landfast ice. In this article we test whether adding tensile strength to sea ice rheologies is able to improve landfast sea ice modeling.
In section 2 the viscous–plastic (VP) rheology introduced by Hibler (1977) and the elastic–viscous–plastic (EVP) rheologies by Hunke and Dukowicz (1997) and by Hunke (2001) are modified to allow for tensile strength. In section 3 the model to test the effect of those changes is described; section 4 discusses the results. We end with our conclusions in section 5.
2. Landfast sea ice rheology
The force balance in the ice pack is given by a two-dimensional momentum equation described further below in Eq. (14). To close the system of equations, a sea ice rheology is needed that allows one to write the internal ice stress σ in terms of the basic variables describing the ice behavior, like sea ice deformation (or strain rates) or ice mechanical properties (e.g., Tremblay and Mysak 1997). In the following the viscous–plastic rheology introduced by Hibler (1979) is rederived to allow for tensile strength (i.e., resistance to divergence). This means that the elliptical yield curve will reach into to the first quadrant of the principal stress diagram (see Fig. 1).
The maximal compressive strength is denoted P as usual, while the maximum tensile stress that the ice can withstand is T. Both numbers are assumed to be positive (i.e., P is positive for compressive strength while T is positive for tensile strengths; see Fig. 1).
The center of the ellipse is thus at (−(P − T )/2, −(P − T )/2). The ratio of long to short semiaxis of the ellipse is e, which normally is set to 2 (Hibler 1979). The semiaxis lengths become 2(P + T )/2 and 2(P + T )/(2e).
The modified equation for the compressive–tensile elliptical yield curve is
In terms of the components of the stress tensor (σ11, σ22, and σ12 = σ21) the same equation looks as follows:
a. Constitutive law with tensile strength
We assume the ice to obey the two-dimensional constitutive law (e.g., Hibler 1979)
where σij is the stress tensor; ε̇ij is the strain rate tensor; ζ and η are bulk and shear viscosities, respectively; p/2 is a pressure term that depends on the local ice properties; and δij is the Kronecker delta function.
Solving Eq. (3) for the strain rates, we get
According to the normal flow rule (Goodier and Hodge 1958), it is assumed that the ice breaks in the direction perpendicular to the yield curve and that the principal axes of strains coincide with the principal axes of stress. Hence, ε̇ij = γ(∂F/∂σij), where F had been defined in Eq. (2) and γ is a proportionality constant. Assuming that the σijs are independent from each other, we take the derivatives of F with respect to the stresses and get
Solving this equation for the σijs we can resubstitute the result into Eq. (2) and find γ to be
where Δ is being introduced for brevity.
Comparing this to the general constitutive law Eq. (4), we see that all the terms line up nicely and thus we can define the bulk and shear viscosities and the pressure term as follows:
Small strain rates can make Δ arbitrarily small, letting the viscosities become infinite. Hibler (1979) regularized this behavior by bounding them by the large values ζmax = 2.5 × 108 g s−1 × P and ηmax = (ζmax/e2). This is equivalent to bounding Δ by a minimal value of Δmin = 2 × 10−9 s−1.
b. Adding elasticity
One way to numerically solve the equations presented so far is to add an elastic term, as was proposed by Hunke and Dukowicz (1997) and further refined in Hunke (2001). The constitutive law Eq. (4) then turns into
where E corresponds to Young’s modulus. While sea ice has a certain elasticity on small scales and such a term makes sense physically (Coon et al. 1974; Colony and Pritchard 1975; Pritchard 1975, 2001; Holland 2007), the values of E used here are purely motivated by its numerical effect and have nothing to do with the elastic property of sea ice. Nonetheless, the modifications assume a new elastic–viscous–plastic rheology of sea ice.
The idea behind adding the elastic term is to find the equilibrium solution (i.e., to have the time derivative go to zero over time). Once the elastic term has approached zero, Eq. (11) becomes the original viscous–plastic constitution law and we recover the desired solution of the viscous–plastic rheology. The Young’s modulus has thus to be chosen in a way that any elastic waves have enough time to die away before the forcings and other variables are updated.
In Hunke and Dukowicz (1997), E is defined as 2E0ρihi(Δx/Δte)2, where E0 is an elasticity parameter that should be chosen between 0 and 1, ρi the density of sea ice, hi the ice thickness, Δx the size of the grid cell, and Δte the time step of the dynamics subcycle. In Hunke (2001), on the other hand, E is defined in terms of a damping time scale for elastic waves Tdamp as E = (ζ /Tdamp) (Hunke 2001, where the parameter is called T ), where ζ is the bulk viscosity of sea ice.
The other difference between the two formulations was that Hunke and Dukowicz (1997) held the viscosities constant over the subcycling to approach as closely as possible the results from Hibler (1979). Hunke (2001) showed that updating the viscosities within the subcycle improved the behavior of the ice model. Our implementation of the two rheologies updates the viscosities within the subcycling. The differences shown thus stem from the differences in the definition of the elasticity parameter E alone.
Adding the elastic term also has the effect that now explicit numerical methods can be used with larger time steps. The viscous–plastic rheology Eq. (3) itself requires much smaller time steps, which make it computationally more expensive.
In areas of slow ice movements, the viscosities can grow very large, forcing the explicit time steps to become very small. The introduction of the elasticity alleviates this situation and allows much bigger time steps. A detailed analysis of this is given in Hunke and Dukowicz (1997).
As suggested by Hibler (1979), we set P to
with the constants as shown in Table 1.
We choose the tensile strength T to be proportional to P. Here P changes greatly with changing ice properties, and hence so does the size of the ellipse. Choosing T proportional to P guarantees that the relative position of the ellipse to the origin remains the same. For example, if T is chosen as constant, the ellipse could move fully into the first quadrant when P becomes small:
Here, kT is set to values in [0, 1]. Setting it to 0 recovers the regular ice rheology, whereas setting it to 1 assumes that ice can resist as much tension as it can compression. Tremblay and Hakakian (2006) estimate values of 0.5 to 0.8 for kT from their analysis of satellite-derived sea ice drift maps.
3. Description of the model
a. Modeling assumptions
Sea ice normally appears in a wide variety of thicknesses, floe shapes, and sizes. To make modeling tractable, the ice is considered to be a two-dimensional continuum (e.g., Pritchard 1975; Hibler 1980; Gray and Morland 1994). We further assume that two variables define the sea ice cover: a is the ice concentration and h is the equivalent ice thickness (Hibler 1979).
The ice concentration a is the fraction of the area in a grid cell covered by ice and thus will go from 0 to 1 (i.e., from 0% to 100%). Values of a below 1 imply that subgrid-scale-sized leads have opened so that 1 − a of the area is open ocean.
Following Hibler (1979), we define the equivalent ice thickness h to be the total volume of ice in a grid cell divided by the area of that grid cell. Thus, h is proportional to the total volume and total mass of ice (as we assume a constant ice density) in a grid cell but has a unit of length. The mass of sea ice in a cell is ρihΔxΔy, where ρi is the density of sea ice and Δx and Δy are the horizontal dimensions of the grid cell. This definition of h is different from the one used by some others, for example by Hunke and Dukowicz (1997).
b. Momentum conservation
The momentum balance of the two-dimensional movement of ice is
The term on the left side describes the change of momentum, with m being the ice mass per unit area and u the ice velocity. The nonlinear advection term has been neglected here.The first term on the right side describes the Coriolis effect, with the Coriolis parameter f and the unit vector normal to the surface k. The atmosphere and ocean stresses are τa and τo, respectively. Finally ∇ · σ is the divergence of the ice stress tensor; this term describes the forces stemming from ice interactions like rafting, ridging, and fracturing.
1) Atmosphere and ocean stresses
Atmosphere and ocean stresses are of the form
where a subscript j would be a for the atmosphere or o for ocean stresses; ρa and ρo are the atmosphere and ocean densities, respectively (assumed to be constant); ca and co are the air and water drag coefficients (see Table 1 for the values used); the wind, ocean current, and ice velocities are ua, uo, and u; and a is the ice concentration (area fraction) and is included here following the suggestions by Gray and Morland (1994) and Connolley et al. (2004). Turning angles are ignored as we use surface wind forcing.
Wind velocities are generally much larger than ice velocities so that we assume ua ≈ ua − u. The equation for the wind stress becomes
Ocean effects are ignored in this study except for a passive drag equivalent to the ice lying on top of an ocean at rest (uo ≡ 0 at all times). The ocean drag term used in this model thus is
c. Conservation laws
From our definitions of the ice concentration a and equivalent ice thickness h, it becomes clear that the following conditions have to be fulfilled at all times:
In our present framework, we ignore thermodynamic processes (i.e., we will assume that we have no freezing or melting). Hence, the total mass of sea ice has to be conserved by the dynamics. Because we assume the ice density to be constant, this is equivalent to assuming ice volume conservation or the conservation of the equivalent ice thickness h:
Additionally, we require that the average ice thickness h/a be conserved except for ridging or rafting. The “average ice thickness” is the thickness of ice averaged over the area covered by ice only; h, on the other hand, is the average over the entire grid cell area. Together with the conservation of h, the conservation of h/a can be reformulated as a conservation of area fraction:
Because the ice concentration a is restricted to [0, 1], it is not truly conserved. Once 1 has been reached under continued convergence, a remains constant while the mass will continue to increase. The average ice thickness h/a is increasing too in such a case; that is, the ice is being crushed and getting thicker, which approximates correctly the behavior of ridging and rafting.
d. The one-dimensional equations
To demonstrate the effect of the modified rheology, we implement a horizontally one-dimensional model (purely zonal flow). We assume υ = 0 and no variation in the y direction (Gray and Morland 1994; Hunke 2001). Most variables retain their meaning; of note is the stress tensor that turns into a single stress variable σ, which corresponds to σ11 of the two-dimensional case (i.e., stress in the x direction applied to a surface pointing into the x direction).
In contrast to the momentum Eq. (14), we ignore the Coriolis force to maintain a purely one-dimensional motion. All we retain are the ice rheology and the forcing terms:
where u is the ice velocity.
There is no shear stress in one dimension and any derivatives with respect to y vanish. The elliptical yield curve collapses to a line. The constitutive law Eq. (3) (i.e., the equation for the stress σ) becomes
where the bulk viscosity ζ = (P + T )/(2Δ), the compressive strength P, and the tensile strength T = kTP have been defined above. From Eq. (6) Δ simplifies considerably to become
Conservation of the equivalent ice thickness h and the ice concentration a become
where 0 ≤ a ≤ 1.
This provides us with the necessary four equations to solve for the four unknowns u, σ, h, and a.
e. Numerical implementation
We solve the four Eqs. (22), (25), (27), and (28) on a staggered grid (one-dimensional Arakawa C grid; see Fig. 2). Ice stress, thickness, and concentration are defined at one set of grid points (the material points), whereas ice velocities are defined in between and to the left and right on the velocity points. Thus, there is one more velocity point than there are material points.
To illustrate the indexing used below, we explain the expression uhik. Except for ρi (density of sea ice), subscripts i to the right indicate the location of a variable (i.e., the grid index). Note that the actual location of a variable with subscript i depends on the variable as we use a staggered grid. Superscripts k to the right indicate the time index. The current (known) time level is k; the next time step is at k + 1. Because the momentum equation is solved on the velocity points, variables defined on the material points such as h have to be interpolated (i.e., averaged) onto the velocity points. Such interpolated values are denoted with a superscripted u on the left side of the variable. Thus, uhik indicates the equivalent ice thickness h at time k interpolated onto the ith node of the velocity grid.
On the left side of our domain we enforce a “closed” boundary condition, simulating the coastline. The velocity along the closed boundary is set to and kept at zero (Dirichlet boundary condition). On the right side there is an “open” boundary that allows ice to leave the domain. This is implemented by enforcing zero derivatives of velocity at the boundaries (Neumann boundary condition) by setting the velocity on the boundary equal to the velocity just inside the domain. Furthermore, the ice strength P is kept at zero on the open boundary.
To close the momentum equations we need to make an assumption about how the internal ice stress depends on other ice variables. In the following subsection the different ice rheologies are explained.
1) Solution of the momentum equation using the viscous–plastic rheology
We use an implicit discretization so we can use larger time steps. Taking first-order differences in space as well as in time (e.g., Durran 1999) leads to a second-order accuracy scheme because of the staggered grid. The momentum Eq. (22) becomes
We use τak = caρauak|uak|uak for brevity as it is a purely external forcing term.
The discretized constitutive law Eq. (25) is
Substituting σ into the discretized momentum equation, we end up with a system of equations of the form
where 𝗔 is a tridiagonal matrix and b contains ice velocities and properties at time k.
To solve, uk+1 = (𝗔k)−1bk is iterated with updating 𝗔 and b in between until the changes in the 2-norm of u (the square root of the sum of the squares of the velocities at each grid point) are less than 10−6. This ensures that the change of a velocity over an iteration is at most 1 μm s−1. In this physical context such a change has become unmeasurable.
2) Solution of the momentum equation using the elastic–viscous–plastic rheology
We use first-order differencing in space as well as time (Euler forward), again leading to a second-order scheme because of the staggered grid.
To find the elasticity parameter E = 2E0ρihi(Δx/Δte)2, as suggested by Hunke and Dukowicz (1997), we can select E0 as well as Δte while the other parameters are given by physics and the modeling results. Following Hunke and Dukowicz (1997) we set E0 = 0.25 (see Table 1). The elastic time step Δte itself is N times smaller than the advective time step Δt. Note that N can be chosen freely and corresponds to the number of iterations that the momentum and rheology equations are iterated. Hunke and Dukowicz (1997) used values from 72 to 400 for N; we use 250. The advective time step Δt is identical to the time step used in the implicit scheme described above. Only every advective time step the ice thickness and concentration are updated.
In Hunke (2001) ,E is defined in terms of a damping time scale for elastic waves Tdamp. For stability, Tdamp has to be chosen between the elastic time step Δte and the advective time step used in the momentum equation Δt. The widely used Los Alamos sea ice model (CICE) uses Δte:Tdamp:Δte = 1:43.2:120 (Hunke and Lipscomb 2008), meaning that Tdamp = 0.36Δt and that N = 120 subcycles are used for the EVP dynamics model. The new version of the Louvain-la-Neuve sea ice model LIM3 (Vancoppenolle et al. 2008) uses Δte:Tdamp:Δte = 1:100:300. We use Δte:Tdamp:Δte = 1:83.3:250, which lies in the same range (Table 1).
Hunke (2001) points out that the updated EVP encounters stability problems in regions of rigid ice (his section 6.1.2). As landfast sea ice is certainly rigid, this has to be addressed. To prevent those instabilities, the following inequality has to be true [Hunke 2001, his Eq. (25)]:
with C being a tuning parameter. Implementations of EVP01 in CICE and LIM3 do not enforce this restriction by default. To investigate its effect on landfast sea ice modeling we used EVP01 in two configurations, once without enforcing Eq. (31) and once with tuning C. For our one-dimensional case, we replaced Δy with Δx and found the optimal C to be 1000 (Table 1).
3) Conservation of equivalent ice thickness and ice concentration
Ice thickness and ice concentration conservation are guaranteed using a first-order upwind scheme. Should the ice concentration have become bigger than 1, it is reset to 1. Both variables are bounded below by 10−6. This is to prevent the ice mass in a grid cell from approaching zero as this would cause numerical problems.
The default model specifications are listed in Table 2. This sets up a strip of landfast ice of 100 km width that is attached to a coast (the closed left boundary) that can leave the domain through the open right boundary should it break loose, although we only analyze results up to break-off. An offshore wind of 10 m s−1 forces the ice away from the coastline. This is probably stronger than what is observed, especially over time scales of weeks, but it serves well to explore the capabilities of models. The tensile strength is set to equal compressive strength, which is rather on the high side. More realistic values seem to find values for the tensile strength around 0.5–0.8 times the compressive strength (Tremblay and Hakakian 2006) for landfast sea ice, which might go close to zero for pack ice, depending on floe size and grid resolution. This choice of setting T = P was made here to conceptually illustrate the effect of adding significant tensile strength.
g. Experimental setups used
We are exploring two different effects: adding tensile strength, on the one hand, and four different numerical formulations, on the other hand. Simulations without tensile strength are called standard, whereas when using the rheology with added tensile strength we call it the tensile setup.
We call the different numerical formulations VP for the viscous–plastic rheology (Hibler 1979), EVP97 for the elastic–viscous–plastic rheology suggested in Hunke and Dukowicz (1997), and EVP01 for the newer formulation of the elastic–viscous–plastic rheology (Hunke 2001). We also look at two different configurations of EVP01. First, we do not enforce the inequality (31) as is done in most implementations, in the EVP01-noC setup. Second, we enforce the inequality with a tuned C in the EVP01-wiC setup. Overall we compare eight setups: VP, EVP97, EVP01-noC, and EVP01-wiC, each as -standard and -tensile.
a. Analytical solutions
We derive some analytical solutions to which our numerical results can later be compared.
1) Maximum width of landfast ice
Assuming we have a sheet of landfast ice of length L (see Fig. 3) at rest, we can solve the momentum Eq. (22) analytically. The ice is at rest everywhere and thus the ocean drag is zero. The equation simplifies to
that is, the stress decreases linearly within the ice. With an ice concentration of 1 within the ice, and an idealized offshore (i.e., positive) wind that is constant over time and space, the wind stress simplifies to
over the landfast ice.
We want to find the maximum width Lmax of landfast ice that can resist the stress exerted by winds of a certain velocity ua. The maximum tensile stress that the ice can sustain is T = P, which will be reached at x = 0:
Using parameters as given in Table 1 and standard values as given in Table 2 (in particular an offshore wind of 10 m s−1, an ice concentration of 1, and an equivalent ice thickness of 1 m.). We find that landfast ice should be able to sustain itself up to a width of Lmax = 212 km. In this calculation, we ignored that the ice starts to move slowly because of ice creep.
2) Ice creep
Ice as modeled by the viscous–plastic rheology undergoes a slow creep when experiencing a stress between the maximum compressive or tensile stress. From the one-dimensional constitutive law Eq. (25) we find that in the viscous regime the strain rate depends on the stress in the following way:
where σ(x) is the stress profile within the landfast ice [see Eq. (34)]. Again choosing standard values, and in particular T = P, the expression simplified considerably (P − T = 0 and P + T = 2P).
The ice velocity at the coastline is 0, so we can find the velocity at any location x in the landfast ice by integrating ∂u/∂x away from the coast:
which is a parabola that goes through the origin and achieves its vertex at the ice edge (see Fig. 4). The modeled creep velocity at the ice edge becomes
Using typical parameter values (see Table 2), the velocity at the ice edge becomes u(100 km) = 1.89 × 10−6 m s−1.
b. Modeling results
1) Effect of adding tensile strength
Figure 5 shows the effect of adding tensile strength in the VP implementation. The images on the left side show the ice thicknesses h of the VP-standard setup, while the right side shows that the ice remains fast in the VP-tensile setup.
2) Comparison of the rheologies
The results from using the VP-tensile setup are very close to the analytical solutions found earlier. Figure 4 shows the velocity profile within the landfast sea ice after one time step. It approximates the analytical solution Eq. (38) to closer than 10−6 m s−1. Note that the momentum equation had been iterated until the solution changed by less than 10−6 m s−1.
The solution from using the three EVP-tensile setups looks quite different (Fig. 6). The ice is able to resist the offshore forcing to varying degrees, but after a few days at the latest it has broken off. In the EVP97-tensile setup, sea ice manages to remain fast for up to 3.5 days. Using the EVP01-noC-tensile, sea ice breaks off almost immediately. Additionally, numerical noise can be seen, which is caused by the instabilities in rigid ice conditions (Hunke 2001). EVP01-wiC-tensile shows some improvement stemming from using the inequality (31) and is able to sustain landfast sea ice for up to a day. Additionally, it can be seen that the numerical noise has been dampened out.
A closer look at the development of ice thickness throughout the landfast ice is given in Fig. 7, where the leftmost column shows the slow thinning of the ice due to the viscous creep as found in the VP-tensile setup. The thinning in the other three columns (EVP-tensile setups) is on a much larger scale. (Note the different scales used on the vertical axis.) In the EVP97-tensile setup (second column) at around 10 km an initially small drop magnifies until it leads to ice break-off after 3 days.
The EVP01-noC-tensile setup clearly suffers from numerical noise. As anticipated, problems appear in this area of rigid ice that lead to a very quick break-off. Enforcing inequality (31) efficiently dampens any elastic waves and improves landfast sea ice modeling to some degree. However, an increased creep leads to accelerated thinning of the ice, which in turn leads to a weakening of the ice strength. This increases the creep rate further, as indicated by the inverse relationship of u(x) with the ice strength P in Eq. (37).
A closer look at the solution from the EVP97-tensile setup shows large velocity fluctuations during the first few time steps as shown in the top of Fig. 8. Those fluctuations do die away over time and the velocity profile approaches the shape of the analytical solution except that the velocities are larger by about a factor of 100, on the order of 10−3 m s−1. These are still very small velocities and probably insignificant in regular geophysical settings, but in this case they change the capability to sustain landfast sea ice drastically.
The initial fluctuations lead to significant initial weakening of the ice floe (see bottom of Fig. 8). Once the ice thickness and concentration (not shown) are diminished, the strength of the ice sheet is significantly reduced, which leads to faster creep rates, as indicated by the inverse relationship of u(x) with the ice strength P in Eq. (37). Increased creep again leads to an accelerated thinning of the ice, thus accelerating the decay. These factors lead to the modeled landfast ice breaking off much earlier than theoretically expected.
The performance of all the EVP-tensile setups can be considerably improved by either increasing N, the number of subcycles of the dynamics model, or by decreasing the size of the advective time step (results not shown). For example, decreasing the advective time step by a factor of 10 also decreases the ice velocities by about a factor of 10 and the ice decays more slowly and remains fast longer. Increasing N by a factor of 4 to 1000 has a similar positive effect on the modeled sea ice. On the other hand, reducing N to 120 (as is used in CICE) further degrades the model performance.
3) Influence of width of landfast sea ice
Finally, we explore the performance of the VP-tensile and the three EVP-tensile setups with respect to different widths of initial landfast sea ice. In particular, we look at how long landfast sea ice of a certain width remains fast under the same offshore wind used above. Figure 9 shows for how many days landfast sea ice of different initial widths remained fast (to an accuracy of 6 h and up to a maximum of 40 days).
The solution of the VP-tensile setup again approaches the theoretical expectation in that ice up to around 212 km remains fast over an extended period of time. The EVP97-tensile setup manages to sustain about a third of the expected width of landfast ice over a significant time duration. Using the EVP01-noC-tensile setup, we find once more that it is not capable of sustaining landfast sea ice. Adding the damping constraint as done in the EVP01-wiC-tensile setup improves matters considerably and small extents of landfast ice manage to remain attached to the shore.
We tested whether adding tensile strength to widely used sea ice rheology formulations enables modeled sea ice to exhibit landfast sea ice features, and in particular to resist movement when subjected to strong offshore winds. We examined the viscous–plastic (VP) rheology by Hibler (1979) as well as elastic–viscous–plastic (EVP) rheologies as described by Hunke and Dukowicz (1997) and modified by Hunke (2001). For this last case, we compared two versions, one in which we do not enforce inequality (31) (EVP01-noC setup) and one in which we did (EVP01-wiC). In each case we compared the tensile setup, where a tensile strength equal to the ice’s compressive strength had been added, to the standard setup without tensile strength.
The solution of the VP-tensile setup does almost as well as our theoretical solution predicted and hence approaches the analytical solution closely. When using the three EVP-tensile setups, on the other hand, the modeled sea ice has a tendency to break off significantly faster. EVP97-tensile manages to hold sea ice of over 50 km fast for several weeks. While this is far from the theoretical expectations, it seems promising enough to test in a more realistic ice-ocean model. EVP01-noC-tensile suffers from a lot of noise that renders it incapable of maintaining fast ice–like features. EVP01-wiC-tensile shows that enforcing inequality (31) as suggested in Hunke (2001) for regions of rigid ice indeed clearly improves the modeled behavior. For application in more realistic models, enforcing this inequality certainly would be necessary. However, based on our results, we conclude that this might be insufficient for the fast ice regions and lead to undesirable consequences elsewhere at the same time.
The performance of all the EVP-tensile setups can be considerably improved by either increasing N, the number of subcycles of the dynamics model, or by decreasing the size of the advective time step. Both changes allow the model solution to converge closer to the analytical solution and increase the time ice remains fast. In realistic applications, increasing N or reducing the size of the time step is computationally expensive and quickly becomes unaffordable. One of the main advantages of EVP over VP, computational efficiency, would be lost.
In all cases the ice always breaks at or very close to the coastline. This is expected from our model because the ice stress is largest there and it is there where the atmospheric pull would first exceed the tensile strength of the ice. In reality, breaking can take place anywhere and often is observed farther away from the coast (Wadhams 1986; Barber and Hanesiak 2004; Mahoney et al. 2007a). It is likely that other effects—for example, grounded ridges, wave action, tides, collisions with pack ice, or thermodynamics—play a role. Depending on the geometry of the coast line (e.g., in narrow strait or small bays), shear strength can also be an important mechanism.
Ultimately, to decide if adding tensile strength to an ice rheology improves the modeling of landfast sea ice, one has to run the modified model including ice thermodynamics on a realistic domain with realistic atmospheric and oceanic forcing and compare the outcome with high-resolution sea ice data. Also, such a framework is needed to more conclusively compare the performance of the various rheologies. The first of such trials by adding tensile strength in a finite element model with high resolution along the coastlines (see description in Lietaer et al. 2008) failed to produce landfast sea ice–like features. This probably resulted from using the EVP01-noC-tensile setup, that is, not enforcing inequality (31). Other approaches such as using EVP01-wiC-tensile are part of ongoing research.
Another open question is the effect of added tensile strength to the ice flow in the ice pack. Miller et al. (2005) found that adding tensile strength by lowering the eccentricity of the ellipse improved the ice thickness distribution across the Arctic. This indicates that adding tensile strength might have a positive effect even away from the coasts. Depending on the model and region, it might also deteriorate modeling results and a combined model might be appropriate—one that treats ice differently in the pack and close to land. Some measurements have indicated that some physical parameters differ fundamentally between pack and fast ice (Prinsenberg et al. 1997) so that such a “phase change” might be needed to properly model landfast sea ice within a conventional pack-ice modeling framework. How such criteria would be formulated is the subject of continuing research.
This work was supported by NSF Office of Polar Programs Grants 07-31682 and 03-37073. We thank the two anonymous reviewers for their helpful comments.
Corresponding author address: Christof König Beatty, Université Catholique de Louvain (ASTR), Chemin du Cyclotron 2, 1348 Louvain-la-Neuve, Belgium. Email: firstname.lastname@example.org