A hybrid ice-mélange model based on particle and continuum methods

Ice mélange, a composite of sea ice and icebergs, can have a major influence on sea-ice–ocean interactions. However, ice mélange has not been represented in climate models because numerically efficient realizations do not exist. This motivates the development of a prototype of a dynamic hybrid ice-mélange model that we present in this paper. In our approach, icebergs are included as particles, while sea ice is treated as a continuum. To derive a joint continuum for ice mélange, we integrate particle properties into the sea-ice continuum. Thus, icebergs are viewed as thick, compact pieces of sea ice. The ice-mélange formulation is based on the viscous–plastic sea-ice rheology, which is currently the most widely used material law for sea ice in climate models. Starting from the continuum mechanical formulation, we modify the rheology such that icebergs are held together by a modified tensile strength in the material law. Due to the particle approach, we do not need highly resolved spatial meshes to represent the typical size of icebergs in ice mélange (< 300 m). Instead, icebergs can be tracked on a subgrid level, while the typical resolution of the sea-ice model can be maintained (≥10 km). This is an appealing property for computational efficiency and for an inclusion within large-scale models. In idealized test cases, we demonstrate that the proposed changes in the material law allow for a realistic representation of icebergs within the viscous–plastic sea-ice rheology. Furthermore, we show that subgrid dynamics, such as polynya formation due to grounded icebergs, can be modelled with the hybrid approach. Overall, this suggested extension of the viscous–plastic sea-ice model is a promising path towards the integration of ice mélange into climate models.

1 Introduction

Fjords with marine-terminating glaciers are commonly found in the polar regions, e.g. around Greenland. These fjord systems can be filled with sea ice into which icebergs calve so that a mixture of sea ice, bergy bits and icebergs is formed: ice mélange. The ice mélange consists of many interacting small icebergs (<300 m) (Dowdeswell et al., 1992; Sulak et al., 2017).

Observations based on field campaigns and remote sensing data indicate that ice mélange affects the glacier–fjord system either by releasing fresh water into the fjord (Enderlin et al., 2018; Moon et al., 2018; Mortensen et al., 2020) or by creating a force at the glacier termini (Cassotto et al., 2015; Bevan et al., 2019; Xie et al., 2019). This force might have enough strength such that it prevents calving events (e.g. Amundson et al., 2010; Krug et al., 2015; Bassis et al., 2021), whereas the release of fresh water through icebergs influences the fjord circulation and melting at the glacier termini (Davison et al., 2020).

Observing ice mélange is difficult, due to the sparsity of remote sensing data and due to the challenges of taking in situ measurements in dense ice conditions. To obtain insights into the potential impact of ice mélange on glacier calving and the underlying ocean circulation, numerical models are therefore necessary. The existing approaches used to include ice mélange into models vary in their complexity. They range from the full description of the ice-mélange dynamics (Robel, 2017; Burton et al., 2018) to the parameterization of specific interactions such as the load which the ice mélange creates at the glacier termini (Schlemm and Levermann, 2021).

There are two distinct approaches in the scientific literature to model the dynamics of ice mélange. In particle methods the ice mélange is expressed using discrete interacting particles, with single particles representing icebergs or sea-ice floes (e.g. Robel, 2017; Burton et al., 2018). As discussed by Vaňková and Holland (2017), implementing a full particle approach into climate models would be extremely challenging due to the enormous numerical costs. In contrast, in the continuum approach, the ice mélange is prescribed as a single continuum (e.g. Pollard et al., 2018; Amundson et al., 2024).

To provide a simple coupling between ice-mélange and sea-ice modules used in climate models, Vaňková and Holland (2017) introduced a continuum ice-mélange model, where sea ice and icebergs build a joint continuum. Icebergs are represented via thick compact pieces of sea ice, which are tracked with a Lagrangian advection using moving meshes. In their model, the icebergs are held together via a large tensile (and shear) strength, which is introduced by a modification of the underlying cavitating fluid sea-ice rheology (Flato and Hibler, 1992). The approach of Vaňková and Holland (2017) requires a high spatial resolution to resolve icebergs in the ice mélange. It is challenging to efficiently solve the nonlinear momentum equation of the underlying sea-ice model with existing solvers (Koldunov et al., 2019; Mehlmann and Richter, 2017b). Currently, efficient solvers for the ice-mélange model of Vaňková and Holland (2017) are missing.

To overcome these difficulties, we introduce a hybrid ice-mélange model. In this approach, the ice mélange is described as a joint continuum consisting of sea ice (continuum) and icebergs (particles). The use of particles in the hybrid approach allows us to track the icebergs on a subgrid level. This has the advantage that icebergs do not need to be explicitly resolved by the spatial mesh. Thus, the typical grid size of several kilometres for a sea-ice model can be used to simulate ice mélange. We derive the momentum equation of the ice mélange by selectively modifying the tensile strength of the sea-ice rheology. This concept is similar to the approach of Vaňková and Holland (2017), but instead of applying the cavitating fluid sea-ice rheology, we consider the viscous–plastic (Hibler, 1979) material law, which has been shown to be more realistic than the cavitating fluid model (Kreyscher et al., 2000).

So far, most climate models have treated sea ice as a viscous–plastic material using the viscous–plastic (Hibler, 1979) or elastic–viscous–plastic (Hunke and Dukowicz, 1997) sea-ice rheology. These rheologies are used in 30 out of the 33 global climate models of the Climate Model Intercomparison Project 5 (CMIP5) (Stroeve et al., 2014). Furthermore, the study of Amundson and Burton (2018) indicates that ice mélange exhibits viscous–plastic deformations. Thus, an inclusion of ice mélange into climate models via a modification of the viscous–plastic material law is a promising approach.

The paper is structured as follows. Sect. 2 presents the ice-mélange model, and Sect. 3 outlines the used numerical discretization. The model is numerically evaluated in Sect. 4. We discuss our results in Sect. 5 and summarize our conclusions in Sect. 6.

2 Viscous–plastic ice-mélange model

Based on the representation of the sea-ice dynamics with a viscous–plastic rheology, we develop a model for ice mélange. As in the underlying sea-ice model (Hibler, 1979), the ice-mélange model consists of three prognostic variables: the ice-mélange's thickness H, its concentration A within a specific grid cell and horizontal velocity v. Ice mélange is considered a joint continuum of sea ice and an iceberg distribution, integrated from a set of iceberg particles {p}. On the continuum level, the icebergs are interpreted as thick and compact pieces of sea ice. Every iceberg is represented by a disc-shaped particle p, which is described by a radius r_p and a height h_p, which can vary between the icebergs. We assume that icebergs are represented by a finite number of small, disc-shaped particles. Thus, the continuum thickness and concentration of ice mélange are described as

$$\begin{array}{}\text{(1)}& {\displaystyle }H(x,y,t)=\left\{\begin{array}{ll}{H}_{\text{sea ice}}(x,y,t)& \text{ if}p(x,y,t)\notin (x,y)\in \mathrm{\Omega },\\ h_p& \text{ if}p(x,y,t)\in (x,y)\in \mathrm{\Omega }.\end{array}\right.\text{(2)}& {\displaystyle }A(x,y,t)=\left\{\begin{array}{ll}{A}_{\text{sea ice}}(x,y,t)& \text{ if}p(x,y,t)\notin (x,y)\in \mathrm{\Omega },\\ \mathrm{1}& \text{ if}p(x,y,t)\in (x,y)\in \mathrm{\Omega }.\end{array}\right.\end{array}$$

Here, Ω is the two-dimensional domain of interest, x and y are the horizontal spatial coordinates, and t is the time. In order to model the ice-mélange velocity, we formulate an expansion of the viscous–plastic rheology that accounts for icebergs. To keep icebergs (thick and concentrated pieces of ice) in the ice-mélange formulation together, we modify the tensile strength of the viscous–plastic law.

In the following, we first generally review the formulation of the governing equations (Sect. 2.1) and the viscous–plastic rheology (Sect. 2.2) before we modify its strength parameterization (Sect. 2.3) to represent icebergs. The iceberg particle interaction and the coupling to the continuum ice-mélange formulation are outlined in Sect. 2.4.

2.1 Momentum and conservation equation

The drift of the ice mélange is described by the two-dimensional momentum equation:

$$\begin{array}{}\text{(3)}& \mathit{\rho }H{\partial }_t\mathit{v}=\mathrm{\nabla }\cdot \mathit{\sigma }+F_{\mathrm{b}},\end{array}$$

where ρ=900 kg m⁻³ is the ice density, and ∇⋅σ describes the divergence of the two-dimensional symmetric stress tensor. The internal stresses are given by the material law described in Sect. 2.3. The remaining terms collected in F_b,

$$\begin{array}{}\text{(4)}& F_{\mathrm{b}}=-\mathit{\rho }Hf\mathit{k}\times \mathit{v}-\mathit{\rho }Hg\mathrm{\nabla }{H}_d+{\mathit{\tau }}_{\mathrm{atm}}-{\mathit{\tau }}_{\mathrm{ocean}}\left(\mathit{v}\right),\end{array}$$

model the body forces acting on the ice mélange: the Coriolis parameter f with upward-pointing unit vector k, the gravity constant g, the surface height H_d, and the atmospheric and oceanic stresses given by τ_atm and τ_ocean. These two drag terms (Coon, 1980) are expressed as

$$\begin{array}{}\text{(5)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{atm}}& {\displaystyle }=C_{\mathrm{atm}}{\mathit{\rho }}_{\mathrm{atm}}{∥{\mathit{v}}_{\mathrm{atm}}∥}_{\mathrm{2}}\left({\mathit{v}}_{\mathrm{atm}}\right),\text{(6)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{ocean}}\left(\mathit{v}\right)& {\displaystyle }=C_{\mathrm{ocean}}{\mathit{\rho }}_{\mathrm{ocean}}{∥\mathit{v}-{\mathit{v}}_{\mathrm{ocean}}∥}_{\mathrm{2}}\left(\mathit{v}-{\mathit{v}}_{\mathrm{ocean}}\right),\end{array}$$

where v_atm describes the wind velocity and v_ocean the ocean current. The corresponding densities are given by ρ_atm=1.3 kg m⁻³ and ρ_ocean=1026 kg m⁻³. The drag coefficients are $C_{\mathrm{atm}}=\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{3}}$ and $C_{\mathrm{ocean}}=\mathrm{5.5}\times {\mathrm{10}}^{-\mathrm{3}}$. Note that ∥⋅∥₂ represents the Euclidean norm. The advection of the sea-ice thickness H_{sea ice} and the concentration A_{sea ice} are calculated as

$$\begin{array}{ll}{\displaystyle }{\partial }_t{H}_{\text{sea ice}}+\mathrm{div}\left(H_{\text{sea ice}}\mathit{v}\right)& {\displaystyle }=\mathrm{0},\\ \text{(7)}& {\displaystyle }{\partial }_t{A}_{\text{sea ice}}+\mathrm{div}\left(A_{\text{sea ice}}\mathit{v}\right)& {\displaystyle }=\mathrm{0},\end{array}$$

with $H_{\text{sea ice}}\in [\mathrm{0},\mathrm{\infty })$ and $A_{\text{sea ice}}\in [\mathrm{0},\mathrm{1}]$. The icebergs are advected based on the continuum ice-mélange velocity v:

$$\begin{array}{}\text{(8)}& {\mathit{x}}_p(t+\mathrm{\Delta }t)={\mathit{x}}_p\left(t\right)+\underset{t}{\overset{t+\mathrm{\Delta }t}{\int }}\mathit{v}\left({\mathit{x}}_p\right(\stackrel{\mathrm{̃}}{t}),\stackrel{\mathrm{̃}}{t})\mathrm{d}\stackrel{\mathrm{̃}}{t},\end{array}$$

where ${\mathit{x}}_p=(x_p,y_p)$ is the position of the centre of particle p.

2.2 Viscous–plastic rheology

Let $\stackrel{\mathrm{̃}}{\mathit{\sigma }}$ be the principal components of the stress tensor σ, given by

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{1}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{11}}+{\mathit{\sigma }}_{\mathrm{22}}}{\mathrm{2}}}+\sqrt{{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{11}}-{\mathit{\sigma }}_{\mathrm{22}}}{\mathrm{2}}}\right)}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{12}}^{\mathrm{2}}},\\ {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{2}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{11}}+{\mathit{\sigma }}_{\mathrm{22}}}{\mathrm{2}}}-\sqrt{{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{11}}-{\mathit{\sigma }}_{\mathrm{22}}}{\mathrm{2}}}\right)}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{12}}^{\mathrm{2}}}.\end{array}\end{array}$$

In the viscous–plastic model (Hibler, 1979), the states of the stress, σ, are described by an elliptic yield curve of the form

$$\begin{array}{}\text{(10)}& F\left({\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{1}},{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{2}}\right)={\left({\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{1}}+{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{2}}+P}{P}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{1}}-{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{2}}}{P}}e\right)}^{\mathrm{2}}-\mathrm{1}=\mathrm{0},\end{array}$$

where e=2 is the minor axis of the ellipse, and P is the ice strength modelled as

$$\begin{array}{}\text{(11)}& P=P^{\star }H\exp (-C(\mathrm{1}-A\left)\right),\end{array}$$

with strength parameter $P^{\star }=\mathrm{27.5}\times {\mathrm{10}}^{\mathrm{3}}$ N m⁻², and C=20. The yield curve is expressed in terms of the principal components of the stress tensor σ:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{1}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{11}}+{\mathit{\sigma }}_{\mathrm{22}}}{\mathrm{2}}}+\sqrt{{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{11}}-{\mathit{\sigma }}_{\mathrm{22}}}{\mathrm{2}}}\right)}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{12}}^{\mathrm{2}}},\\ {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{2}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{11}}+{\mathit{\sigma }}_{\mathrm{22}}}{\mathrm{2}}}-\sqrt{{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{11}}-{\mathit{\sigma }}_{\mathrm{22}}}{\mathrm{2}}}\right)}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{12}}^{\mathrm{2}}}.\end{array}\end{array}$$

A visualization of the yield curve is given in Fig. 1 (blue line). As there are no stress states on the yield curve with $({\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{1}},{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{2}})>\mathrm{0}$, sea ice has no tensile strength (no resistance to divergence) (Leppäranta, 2011). The stress states are related to the strain rates,

$$\begin{array}{}\text{(13)}& {\dot{\mathit{ϵ}}}_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left\{{\partial }_{{\mathbf{x}}_j}{\mathit{v}}_i+{\partial }_{{\mathbf{x}}_i}{\mathit{v}}_j\right\},\end{array}$$

by the constitutive law (Hibler, 1979):

$$\begin{array}{}\text{(14)}& {\mathit{\sigma }}_{ij}=\mathrm{2}\mathit{\eta }{\dot{\mathit{ϵ}}}_{ij}+(\mathit{\zeta }-\mathit{\eta })\left({\dot{\mathit{ϵ}}}_{\mathrm{11}}+{\dot{\mathit{ϵ}}}_{\mathrm{22}}\right){\mathit{\delta }}_{ij}-{\displaystyle \frac{P}{\mathrm{2}}}{\mathit{\delta }}_{ij},\end{array}$$

where δ_ij is the Kronecker symbol. The nonlinear shear ζ and bulk viscosity η are chosen as

$$\begin{array}{}\text{(15)}& \mathit{\eta }=e^{-\mathrm{2}}\mathit{\zeta },\mathit{\zeta }={\displaystyle \frac{P}{\mathrm{2}\mathrm{\Delta }\left(\dot{\mathit{ϵ}}\right)}}.\end{array}$$

To guarantee a smooth transition between the viscous and the plastic regime, we follow Kreyscher et al. (2000) and choose

$$\begin{array}{}\text{(16)}& \mathrm{\Delta }\left(\dot{\mathit{ϵ}}\right)=\sqrt{{\mathrm{\Delta }}_P(\dot{\mathit{ϵ}}{)}^{\mathrm{2}}+{\mathrm{\Delta }}_{\text{min}}^{\mathrm{2}}}.\end{array}$$

In the case of the plastic regime, ${\mathrm{\Delta }}_P\left(\dot{\mathit{ϵ}}\right)$ is defined as

$$\begin{array}{}\text{(17)}& {\mathrm{\Delta }}_P\left(\dot{\mathit{ϵ}}\right)=\sqrt{\left({\dot{\mathit{ϵ}}}_{\mathrm{11}}^{\mathrm{2}}+{\dot{\mathit{ϵ}}}_{\mathrm{22}}^{\mathrm{2}}\right)\left(\mathrm{1}+e^{-\mathrm{2}}\right)+\mathrm{4}{e}^{-\mathrm{2}}{\dot{\mathit{ϵ}}}_{\mathrm{12}}^{\mathrm{2}}+\mathrm{2}{\dot{\mathit{ϵ}}}_{\mathrm{11}}{\dot{\mathit{ϵ}}}_{\mathrm{22}}\left(\mathrm{1}-e^{-\mathrm{2}}\right)}.\end{array}$$

The viscous regime is given as

$$\begin{array}{}\text{(18)}& {\mathrm{\Delta }}_{\text{min}}\left(\dot{\mathit{ϵ}}\right)=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{9}}.\end{array}$$

2.3 Strength parameterization

The absence of tensile strength in the original model is apparent from the fact that the yield curve of the viscous–plastic rheology does not contain combinations of $({\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{1}},{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{2}})>\mathrm{0}$ (see the blue curve in Fig. 1). A tensile strength has been introduced into this model, e.g. by König and Holland (2010) to model landfast sea ice. Similar to Vaňková and Holland (2017), which added a tensile strength into a cavitating fluid model in order to prescribe ice mélange, we introduce a tensile strength into the standard viscous–plastic sea-ice rheology to model icebergs. This tensile strength leads to a resistance to divergence in the presence of icebergs.

Figure 1 Two-dimensional yield curve in principal stress space without (blue) and with (red) tensile strength. Pure divergence d is at the origin of the graph; pure convergence is obtained at c. ${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{1}}$ and ${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{2}}$ are the principal components of the stress tensor σ (see Eq. 14).

By including the tensile strength, the elliptic yield curve is shifted into the first quadrant (red curve in Fig. 1). The new centre of the ellipse is given by $(-\frac{P-T}{\mathrm{2}},-\frac{P-T}{\mathrm{2}})$, with the maximum tensile strength T. Both P and T are positive numbers. Thus, the modified elliptic yield curve is given by

$$\begin{array}{}\text{(19)}& \begin{array}{rl}F\left({\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{1}},{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{2}}\right)& ={\left({\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{1}}+{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{2}}+P-T}{P+T}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{1}}-{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{2}}}{P+T}}e\right)}^{\mathrm{2}}\\ & -\mathrm{1}=\mathrm{0},\end{array}\end{array}$$

with the elliptic ratio e. For T=0 (no tensile strength), the elliptic yield curve is equivalent to the ellipse of the viscous–plastic sea-ice rheology (see Eq. 10). It is assumed that ∂F and the strain rates are perpendicular to the surface of the yield curve (Leppäranta, 2011). Thus, the relation between the stress tensor and the strain rates for the shifted yield curve is derived by a normal flow rule:

$$\begin{array}{}\text{(20)}& {\dot{\mathit{ϵ}}}_{ij}=\mathit{\gamma }{\displaystyle \frac{\partial F({\mathit{\sigma }}_{\mathrm{11}},{\mathit{\sigma }}_{\mathrm{22}},{\mathit{\sigma }}_{\mathrm{12}},{\mathit{\sigma }}_{\mathrm{21}})}{\partial {\mathit{\sigma }}_{ij}}},\end{array}$$

with γ>0. This leads to the modified rheology

$$\begin{array}{}\text{(21)}& {\mathit{\sigma }}_{ij}=\mathrm{2}\mathit{\eta }{\dot{\mathit{ϵ}}}_{ij}+(\mathit{\zeta }-\mathit{\eta })\left({\dot{\mathit{ϵ}}}_{\mathrm{11}}+{\dot{\mathit{ϵ}}}_{\mathrm{22}}\right){\mathit{\delta }}_{ij}-{\displaystyle \frac{P-T}{\mathrm{2}}}{\mathit{\delta }}_{ij},\end{array}$$

with the bulk and shear viscosities as

$$\begin{array}{}\text{(22)}& \mathit{\zeta }={\displaystyle \frac{P+T}{\mathrm{2}\mathrm{\Delta }}},\mathit{\eta }={\displaystyle \frac{\mathit{\zeta }}{e^{\mathrm{2}}}}={\displaystyle \frac{P+T}{\mathrm{2}\mathrm{\Delta }{e}^{\mathrm{2}}}}.\end{array}$$

Following König and Holland (2010), we define the tensile strength T relative to the compressive strength:

$$\begin{array}{}\text{(23)}& T=P^{\star }H\mathrm{\Phi }.\end{array}$$

The indicator function Φ is given as

$$\begin{array}{}\text{(24)}& \mathrm{\Phi }(x,y)=\left\{\begin{array}{l}\mathrm{1}\text{ if }(x-x_p{)}^{\mathrm{2}}+(y-y_p{)}^{\mathrm{2}}<r^{\mathrm{2}},\\ \mathrm{0}\text{ if }(x-x_p{)}^{\mathrm{2}}+(y-y_p{)}^{\mathrm{2}}\ge r^{\mathrm{2}}.\end{array}\right.\end{array}$$

2.4 Iceberg interaction

An interaction of two distinct particles p_i and p_j is modelled by a hard disc model (Herman, 2011) if the particles overlap:

$$\begin{array}{}\text{(25)}& \left|\right|{\mathit{x}}_{p_i}\left(t\right)-{\mathit{x}}_{p_j}\left(t\right)\left|\right|\le r_i+r_j,\end{array}$$

where r_i and r_j are the radii, and ${\mathit{x}}_{p_i}\left(t\right)=(x_{p_i},y_{p_i})$ and ${\mathit{x}}_{p_j}\left(t\right)=(x_{p_j},y_{p_j})$ are the positions of the interacting particles p_i and p_j, respectively. The position of the overlapping particles is corrected by assuming an inelastic collision (Herman, 2011). For this, we use the last particle position at which the particles had not collided and update the location of the particles based on the ice-mélange velocity corrected for the collision. The latter is calculated as follows:

$$\begin{array}{}\text{(26)}& {\stackrel{\mathrm{̃}}{\mathit{v}}}_i={\mathit{v}}_i-{\displaystyle \frac{{\mathit{\alpha }}_{ij}}{m_i}},{\stackrel{\mathrm{̃}}{\mathit{v}}}_j={\mathit{v}}_j+{\displaystyle \frac{{\mathit{\alpha }}_{ij}}{m_j}},\end{array}$$

where ${\mathit{v}}_i=\mathit{v}(x_{p_i},y_{p_i},t)$ and ${\mathit{v}}_j=\mathit{v}(x_{p_j},y_{p_j},t)$ are the current velocities of the ice mélange, and m_i and m_j are the masses of particle p_i and p_j, respectively. The coefficient α_ij is given by

$$\begin{array}{}\text{(27)}& {\mathit{\alpha }}_{ij}={\displaystyle \frac{m_i{m}_j}{m_i+m_j}}(\mathrm{1}+\mathit{ϵ})\cdot n_{ij}({\mathit{v}}_i-{\mathit{v}}_j)n_{ij},\end{array}$$

with the relative unit position $n_{ij}=\frac{{\mathit{x}}_{p_i}-{\mathit{x}}_{p_j}}{\Vert {\mathit{x}}_{p_i}-{\mathit{x}}_{p_j}\Vert }$. The coefficient of restitution ϵ is set to 0.9, as this value has been used in the past for sea ice (Shen et al., 1987). Since we assume that icebergs are thick pieces of sea ice, this choice is appropriate.

3 Numerical discretization

The ice-mélange model is implemented in the open-source academic software library Gascoigne (Braack et al., 2021), which uses quadrilateral grids. On the mesh, the velocity unknowns are placed at the vertices, whereas the tracers are staggered at the cell centres. This placement corresponds to an A-grid- and a B-grid-type staggering for the velocity and tracers, respectively. The velocity is approximated in space with piecewise linear finite elements, whereas the tracers are discretized as a piecewise constant per cell.

For the time discretization of the ice-mélange model, we split the coupled system of equations in time. First, we approximate the solution of the momentum equation (Eq. 3). Then, the solution of the transport (Eqs. 7 and 8) with the updated velocity is computed. This choice of the implicit Euler method is motivated by the fact that an explicit discretization of the viscous–plastic sea-ice model requires a time step of 1 s on a grid with a size of 100 km × 100 km (Ip et al., 1991). We expect similar constraints for the ice-mélange model, because its rheology is based on the viscous–plastic sea-ice model.

For our choice of an implicit temporal discretization, a nonlinear system of differential equations needs to be solved in every time step. We use a modified Newton method for this solution as it shows improved convergence compared to a standard Newton method and Picard solver (Mehlmann and Richter, 2017a).

3.1 Coupling between particle and continuum method

In order to derive the continuum thickness and concentration of the ice mélange, the icebergs in the form of particles need to be numerically coupled into the continuum sea-ice formulation. We realize this by calculating a continuum iceberg thickness and iceberg concentration in each cell K:

$$\begin{array}{}\text{(28)}& {\displaystyle }{A}_{\text{iceberg}}{|}_K=& {\displaystyle }\sum _{p\in K}{\displaystyle \frac{a_p}{\left|K\right|}},\text{(29)}& {\displaystyle }{H}_{\text{iceberg}}{|}_K=& {\displaystyle }\sum _{p\in K}{\displaystyle \frac{h_p{a}_p}{\left|K\right|}},\end{array}$$

with particle area $a_p=\mathit{\pi }{r}^{\mathrm{2}}{|}_K$ and area of a grid cell $\left|K\right|$. We use the iceberg concentration to determine the tensile strength in the presence of icebergs. Thus, the discretized version of the tensile strength is given by

$$\begin{array}{}\text{(30)}& T=\left\{\begin{array}{ll}\mathrm{0}& \text{ if }{A}_{\text{iceberg}}{|}_K<\frac{\mathit{\pi }(\mathrm{0.5}\sqrt{\left|K\right|}{)}^{\mathrm{2}}}{\left|K\right|},\\ P^{\star }H{c}_{\text{tensile}}{A}_{\text{sea ice}}& \text{ else. }\end{array}\right.\end{array}$$

In the discretized version, we account for the sea-ice concentration weighted by a constant c_tensile. For the sake of simplicity, we choose c_tensile=1. The inclusion of the sea-ice concentration allows us, for example, to model no tensile strength between icebergs if there is no sea ice present.

The threshold is selected such that the tensile strength becomes active as soon as an area of the grid cell is filled with icebergs that cover an area at least as large as a disc-shaped iceberg with radius of $\frac{\sqrt{\left|K\right|}}{\mathrm{2}}$. This choice is discussed in Sect. 5.

The ice-mélange concentration and thickness in each grid cell are given by

$$\begin{array}{}\text{(31)}& {\displaystyle }A{|}_K& {\displaystyle }=min(A_{\text{iceberg}}{|}_K+A_{\text{sea ice}}{|}_K,\mathrm{1}),\text{(32)}& {\displaystyle }H{|}_K& {\displaystyle }=H_{\text{iceberg}}{|}_K+H_{\text{sea ice}}{|}_K.\end{array}$$

In the presence of icebergs, sea ice is more compressed and thicker compared to areas without icebergs. The effective sea-ice thickness and sea-ice concentration are given by ${\stackrel{\mathrm{̃}}{A}}_{\text{sea ice}}:=min(A_{\text{sea ice}}/(\mathrm{1}-A_{\text{iceberg}}),\mathrm{1})$ and ${\stackrel{\mathrm{̃}}{H}}_{\text{sea ice}}:=H_{\text{sea ice}}/(\mathrm{1}-A_{\text{iceberg}})$, respectively. In the context of this ice-mélange formulation, we assume that icebergs are represented by a finite number of small particles with $r_p\le \frac{\sqrt{K}}{\mathrm{2}}$.

We summarize the time discretization of the ice-mélange dynamics in Algorithm 1. The time loop starts with the calculation of the iceberg distribution (Eq. 28) and the computation of the ice-mélange tracers (Eqs. 31, 32), which corresponds to step 1 and step 2 of Algorithm 1. Then, the updated ice-mélange tracers are coupled to the momentum equation to solve for the ice-mélange velocity (step 3 of Algorithm 1). To calculate the advection of ice mélange in Eqs. (7) and (8), we separately transport the continuum sea-ice tracers and the iceberg particles (step 4 and step 5 of Algorithm 1, respectively). Sea ice in Eq. (7) is advected via an upwind scheme, while the particles are transported according to Eq. (8) in a substepping procedure. With this approach, each particle is advected with the corresponding ice-mélange velocity. The latter is given by evaluating the piecewise linear finite element interpolation at the particle location. During the substepping procedure, the icebergs are checked for collisions with other icebergs or with the boundary of the domain and then replaced accordingly to Eq. (26).

Algorithm 1 Partitioned time-stepping loop.

Let $I=[\mathrm{0},T]$ be the time span of interest, and let $\mathit{v}\left(t_{\mathrm{0}}\right),\mathit{\left\{}{\mathit{x}}_p\right(t_{\mathrm{0}}\left)\mathit{\right\}},A_{\text{sea ice}}\left(t_{\mathrm{0}}\right)$ and H_{sea ice}(t₀) be the initial solutions of the ice-mélange velocity, the initial position of each particle, the sea-ice concentration and sea-ice thickness at time t₀=0, respectively. The time period is discretized into equidistant steps: $\mathrm{0}=t_{\mathrm{0}}<t_{\mathrm{1}},\mathrm{\dots },<t_N=T.$ The following time iteration is performed for $n=\mathrm{1},\mathrm{2},\mathrm{\dots },N$.

1. Derive an average iceberg concentration A_iceberg(t_n) and iceberg thickness H_iceberg(t_n) (Eq. 28) based on the set of corresponding particle positions {x_p(t_n−1)}.

2. Calculate the ice-mélange concentration A(t_n) (Eq. 31) and ice-mélange thickness H(t_n) (Eq. 32) based on the concentration and thickness of the continuum sea-ice and iceberg distribution, $A_{\text{sea ice}}\left(t_{n-\mathrm{1}}\right),H_{\text{sea ice}}\left(t_{n-\mathrm{1}}\right),A_{\text{iceberg}}\left(t_{n-\mathrm{1}}\right)$ and H_iceberg(t_n−1), respectively.

3. Solve the momentum equation (Eq. 3) based on A(t_n) and H(t_n):

   $$\mathit{v}\left(t_{n-\mathrm{1}}\right)\to \mathit{v}\left(t_n\right).$$

4. Solve the advection Eq. (7) based on the velocity v(t_n):

   $$A_{\text{sea ice}}\left(t_{n-\mathrm{1}}\right)\to A_{\text{sea ice}}\left(t_n\right).$$

5. Based on a partitioning into equidistant substeps $t_{n-\mathrm{1}}=t_{m-\mathrm{1}}<\mathrm{\dots }<t_m=t_n$, calculate for $m=\mathrm{1},\mathrm{2},\mathrm{\dots },M$ the particle position (Eq. 8) based on v(t_n) and account for particle interaction in the case of collision between particles (Eq. 26):

   $${\mathit{x}}_p\left(t_{m-\mathrm{1}}\right)\to {\mathit{x}}_p\left(t_m\right).$$

4 Numerical validation

The proposed hybrid ice-mélange model is tested in six idealized test cases. The first three test cases (Sect. 4.1) highlight the need for introducing a tensile strength for icebergs into the viscous–plastic sea-ice rheology. Based on the verification of the modified tensile strength, we numerically analyse the combination of particle and continuum methods to represent ice mélange in the last three test cases (Sect. 4.2).

4.1 Tensile strength

The first three test cases are designed such that the behaviour of the ice mélange is tested under compressive (Sect. 4.1.1), tensile (Sect. 4.1.2) and shear (Sect. 4.1.3) forces. For the sake of simplicity, we use v_ocean=0 m s⁻¹. In order to reduce the complexity of the analysis of the modified tensile strength, we neglect the particle coupling in these first three test cases. Instead of simulating the iceberg motion via particles and integrating it into the sea-ice continuum, we track the icebergs (thick and compact pieces of sea ice) in the sea ice via an indicator function ϕ that is transported in time:

$$\begin{array}{}\text{(33)}& {\partial }_t\mathit{\varphi }+\mathrm{div}\left(\mathit{\varphi }\mathit{v}\right)=\mathrm{0},\mathit{\varphi }(x,y,\mathrm{0})={\mathit{\delta }}_{xy},\end{array}$$

where δ_xy is the Kronecker symbol that equals 1 in the presence of icebergs. According to the volume-in-fluid method (Hirt and Nichols, 1982), the indicator function Φ of Eq. (23) is modified to

$$\begin{array}{}\text{(34)}& \mathrm{\Phi }=\left\{\begin{array}{ll}\mathrm{0}\text{ if }\mathit{\varphi }\le c& \text{for sea ice,}\\ \mathrm{1}\text{ if }\mathit{\varphi }>c& \text{for icebergs.}\end{array}\right.\end{array}$$

We have chosen c=0.3 by experimental tuning.

4.1.1 Iceberg pushed against a wall

The first test case is similar to the one used by Vaňková and Holland (2017). The domain is given by an area of the size 5 km × 5 km; see left plot in Fig. 2. All boundaries except the exit to the ocean (right boundary) use Dirichlet boundary conditions (v=0). The upper and lower boundaries represent coastlines, while the boundary on the left represents the glacier terminus. We place a 1 km × 2 km large iceberg in front of the glacier wall. A grid cell size of 110 m is used. The iceberg is a 10 m thick compact block of sea ice with a concentration of 1. The left half (<2.5 km) of the domain is filled with 0.1 m thick sea ice with a concentration of 0.1. The iceberg is pushed against the glacier terminus by a constant wind v_atm=20 m s⁻¹.

Figure 2 Visualization of the ice-mélange thickness for a test case in which an iceberg (thick red contour line) is pushed against the glacier terminus (left boundary) by a constant wind (v_atm=20 m s⁻¹). The panels denote the results for (a) the standard viscous–plastic (VP) rheology and (b) for the viscous–plastic rheology with tensile strength. Both snapshots show the simulation result after 3 h.

The simulation is run for 3 h. When comparing the standard viscous–plastic rheology and the modified material law (Fig. 2), we find that for the standard viscous–plastic rheology the thick and compact piece of sea ice, which models the iceberg, deforms towards the glacier and accumulates in front of the glacier terminus. This is visible by the thick red line shown in Fig. 2a, which visualizes the iceberg's contour. The sea ice along the glacier terminus (outside the contour line) piles up. With the proposed modification, the iceberg is able to keep its rectangular shape throughout the simulation (see Fig. 2b). This test case shows that in our modified viscous–plastic rheology, the iceberg can withstand the compression, initiated by the wind forcing.

4.1.2 Iceberg under a divergent wind field

We consider a 5 km × 5 km large domain with a divergent wind field:

$$\begin{array}{}\text{(35)}& {\mathit{v}}_{\text{atm}}=\left\{\begin{array}{ll}\mathrm{15}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}& x>\mathrm{2.5}\mathrm{km},\\ -\mathrm{15}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}& x<\mathrm{2.5}\mathrm{km}.\end{array}\right.\end{array}$$

The 10 m thick iceberg with a concentration of 1 is placed in the middle of the domain. The whole domain is filled with sea ice whose concentration is 0.1 and whose thickness is 0.1 m. Using the standard viscous–plastic rheology, the iceberg is torn apart after 3 h (see Fig. 3a). With additional tensile strength, the iceberg keeps its form (see Fig. 3b). The surrounding sea ice is still transported in the wind direction. This test case shows, in particular, that the iceberg maintains its shape under diverging wind conditions at the location of the iceberg.

Figure 3 The iceberg (thick red contour line) is placed under a divergent wind field. The panels denote the results for (a) the standard viscous–plastic (VP) rheology and (b) for the viscous–plastic rheology with tensile strength. Both snapshots show the ice-mélange thickness after 3 h.

4.1.3 Iceberg under shear

In the third test case, the same initial setup as in Sect. 4.1.2 is used. But instead of using a divergent wind field, we apply a shearing wind field:

$$\begin{array}{}\text{(36)}& {\mathit{v}}_{\text{atm}}=\left\{\begin{array}{ll}\mathrm{15}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}& y>\mathrm{2.5}\mathrm{km},\\ -\mathrm{15}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}& y<\mathrm{2.5}\mathrm{km}.\end{array}\right.\end{array}$$

The domain of size 5 km × 5 km is filled with 0.1 m thick sea ice with a concentration of 0.1.

Figure 4 The iceberg is placed under a shearing wind field (${\mathit{v}}_{\text{atm}}=\pm \mathrm{15}$ m s⁻¹). Panel (a) visualizes the result for the standard viscous–plastic (VP) rheology and (b) for the viscous–plastic rheology with tensile strength. Both snapshots show the ice-mélange thickness after 3 h. The thick red line indicates the iceberg contour.

In both cases, the iceberg rotates clockwise (see Fig. 4) as expected. Using the standard viscous–plastic rheology, the iceberg deforms. Fig. 4a shows that the iceberg contour is slightly S-shaped, and parts from the iceberg already detach. This is in contrast to the behaviour of the iceberg using the additional tensile strength (see Fig. 4b). Here, the iceberg contour stays rectangular, and the iceberg rotates as one sea-ice block.

4.2 A hybrid ice-mélange representation

The advantage of using a particle method on the joint continuum of sea ice and icebergs is given by the fact that the icebergs and their interactions can be modelled on a subgrid-scale. We consider a domain of size 512 km × 512 km covered with a quadrilateral mesh of size 16 km. In this setup, icebergs are represented by particles with a radius of 125 m. Using this configuration, we study an iceberg–iceberg interaction under shear forcing (Sect. 4.2.1), iceberg separation under divergent forcing (Sect. 4.2.2) and the formation of a polynya due to subgrid iceberg grounding (Sect. 4.2.3). All test cases use a time step size of 2000 s.

4.2.1 Iceberg–iceberg interaction

An explicit iceberg–iceberg collision is forced in order to test the behaviour of icebergs under contact. Two icebergs with a height of 20 m are placed into a 2 m thick sea-ice field with a concentration of 0.7. The icebergs approach each other due to the following wind field:

$$\begin{array}{}\text{(37)}& {\mathit{v}}_{\text{atm}}=\left\{\begin{array}{ll}\mathrm{10}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}& y<\mathrm{256}\mathrm{km},\\ -\mathrm{10}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}& y>\mathrm{256}\mathrm{km}.\end{array}\right.\end{array}$$

The velocity field of the ice mélange after the first time step is shown in Fig. 5a. At first, the icebergs approach each other until they collide (Fig. 5b). Due to collision, iceberg 1 moves further up and iceberg 2 further down. After the collision event, both icebergs drift past each other and separate again; see Fig. 5c.

Figure 5 Visualization of the first component, v_x, of the ice-mélange velocity for a test case in which two icebergs are transported towards each other by a constant wind (${\mathit{v}}_{\text{atm}}=\pm \mathrm{10}$ m s⁻¹). Panel (a) shows the whole domain with the initial velocity field. The last two panels show a closeup of the iceberg interaction after (b) 100 time steps and after (c) 200 time steps with a time step size of 2000 s.

4.2.2 Iceberg field under diverging winds

We analyse the response of an iceberg field to diverging winds and consider the following wind field:

$$\begin{array}{}\text{(38)}& {\mathit{v}}_{\text{atm}}=\left\{\begin{array}{ll}-\mathrm{20}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}& x<\mathrm{256}\mathrm{km},\\ \mathrm{20}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}& x>\mathrm{256}\mathrm{km}.\end{array}\right.\end{array}$$

The icebergs used in the setup are 20 m thick. The sea ice, if present, has a concentration of 0.7 and a thickness of 2 m. To highlight the influence of the tensile strength (Eq. 30) on the motion of the iceberg field, we compare different setups in Fig. 6.

Figure 6 Iceberg field under a diverging wind field. The icebergs have a height of 20 m. In panels (a) and (b), sea ice is 2 m thick with a concentration of 0.7, while in panel (c) the icebergs have no surrounding sea ice. All three panels show the simulation after 300 time steps.

In Fig. 6a and b, sea ice is presented between and around the iceberg field. In contrast to Fig. 6b, the tensile strength parameterization is not active in Fig. 6a. Thus, the surrounding sea ice and the icebergs disperse in Fig. 6a. We use the same setup in Fig. 6b but with activated tensile strength. The surrounding sea ice moves apart, according to the wind field, but the iceberg particles stay in their initial grid cell. Without sea ice between and around the icebergs, the tensile strength in Eq. (30) equals zero. Therefore, the icebergs disperse in Fig. 6c. The setup shows that the modified rheology allows for icebergs to disperse if no sea ice is present. At the same time, we see that the tensile strength is necessary to prevent the field from moving in the case of an iceberg field with sea ice in between. The amount of active tensile strength can be controlled by the parameter c_tensile (Eq. 30), and this is further discussed in Sect. 5.

4.2.3 Iceberg grounding

Icebergs can impact sea ice in different ways, e.g. mechanically by colliding and breaking up the sea-ice cover, by creating openings, and by altering the structure of the sea-ice cover. To simulate such an interaction, we simulate iceberg grounding in the final test case. Those grounding events occur in shallow waters and have profound implications for sea-ice dynamics. As icebergs come into contact with the seafloor, they become immobilized, transforming into obstacles that influence the surrounding sea ice. This affects the natural flow and movement of sea ice and thereby the local circulation patterns and the distribution of sea ice. In this test case, we simulate iceberg grounding and analyse the resulting dynamics of the ice mélange with respect to the formation of a polynya. The domain consists of a 2 m thick sea-ice layer with a concentration of 0.7. Three grid cells, each of the size of 16 km × 16 km, are filled with 4096 icebergs per cell (see Fig. 7). Each iceberg has a radius of 125 m and is 20 m thick. We use such a large number of icebergs to show that we can work with many icebergs per cell that are much smaller than the cell size. These icebergs represent the effects of multiple grounded icebergs in the area with different sizes.

The icebergs in the two lower grid cells closer to the boundary are marked as grounded (v_iceberg=0 m s⁻¹). As forcing, we use an ocean current of v_ocean=0.2 m s⁻¹ and neglect any atmospheric forcing (v_atm=0 m s⁻¹).

Figure 7 Panel (a) shows the initial ice-mélange concentration of the grounded iceberg setup (the two lower grid cells). Panel (b) presents a closeup of the upper grid cell so that individual iceberg particles are visible.

Figure 8 Closeup of the ice-mélange velocity after five time steps of the grounded iceberg setup. The area of low velocity in the centre of the domain indicates the presence of the icebergs. The ice mélange is forced by a constant ocean current (v_ocean=0.2 m s⁻¹).

Figure 9 Three grid cells are each filled with 4096 icebergs. The lower two grid cells consist of only grounded icebergs. The snapshots in (a) and (b) show the ice-mélange concentration after 600 time steps with an ocean forcing of v_ocean=0.2 m s⁻¹. Panel (c) shows a closeup of the ungrounded icebergs in (b).

Figure 9a shows the grounding event using the standard viscous–plastic rheology. The ungrounded icebergs in the upper of the three cells are transported as one block of icebergs towards the right boundary. The square shape is slightly deformed and rotated. The surrounding sea ice accumulates in the domain's right half and builds a straight ice edge. This is a nonphysical behaviour as grounded icebergs should cause a pile-up of sea ice in front of the icebergs. Due to the missing tensile strength in this configuration, the ice-mélange velocity is not zero in the dense iceberg field; see Fig. 8a. Thus, the sea ice flows through the icebergs.

This is in contrast to the results conducted with the modified rheology (see Fig. 9b). Here, the additional tensile strength leads to nearly zero velocity in the dense iceberg field (Fig. 8b) and prevents the sea ice from passing through the icebergs. Therefore, sea ice accumulates in front of the icebergs. Compared to the setup with the standard viscous–plastic rheology, the icebergs move only slightly to the right in the configuration with the modified tensile strength. This is visible in the closeup shown in Fig. 9c. The active tensile strength in the upper of the three cells leads to a reduction in the ice-mélange velocity, which results in a smaller displacement of the iceberg particles.

5 Discussion

The results from Sect. 4 show that the combined approach of subgrid iceberg particle dynamics and a continuum formulation is able to simulate ice-mélange dynamics with respect to different wind and ocean forcings. The usage of particle icebergs in the setup allows for a representation of ice mélange on coarse-resolution horizontal meshes, which have the same resolution as the meshes applied for simulating large-scale sea-ice dynamics in climate models.

In line with prior findings, the test cases presented in Sect. 4 demonstrate that a modification of the original sea-ice rheology is crucial to represent icebergs in the viscous–plastic model. In areas with high iceberg coverage, the ice mélange behaves almost like a rigid body due to the modification of the strength parameter.

The ice-mélange model is developed from the sea-ice perspective with the aim to include the dynamical effect of small icebergs on the evolution of the sea-ice dynamics. Therefore, the prescription of the iceberg dynamics in the absence of sea ice is very limited. Icebergs move either due to collision or with an averaged ice-mélange velocity calculated from the ice-mélange momentum equation. One perspective to allow for a more complex motion of icebergs, especially in the absence of sea ice, is to use a particle model with higher fidelity to represent the iceberg motion, e.g. the approach used in Robel (2017).

In order to ensure numerical efficiency, we have represented icebergs on the particle level as round discs. This simplification of the icebergs' geometry effects the simulated iceberg interactions. The use of geometric objects with other shapes can lead to a motion with different directions after the collision. But calculating the collision of more complex geometric objects such as polygons is numerically more expensive compared to the usage of disc-shaped particles (Damsgaard et al., 2021). In addition, the representation via discs requires a uniform iceberg thickness, which may lead to a coarse approximation of the forces in the ice mélange. Since we model the iceberg interaction on a subgrid scale and integrate the icebergs into the large-scale sea-ice model, the impact of these simplifications is of secondary importance.

In the test cases in which we consider iceberg particles, the relative speed between ice mélange and iceberg particles is almost zero. For the sake of simplicity, we neglected the feedback from the modified iceberg velocity to the ice-mélange velocity. We plan to adjust the ice-mélange velocity to account for this feedback. One possibility is to derive an average velocity of the icebergs per cell and include a drag term in the ice-mélange momentum equation that accounts for the difference in the iceberg and the ice-mélange velocity.

Concerning the response of the ice mélange to subgrid iceberg grounding, we note that polynyas, which can be simulated in the ice mélange, cannot be smaller than the size of a grid cell, as in any standard sea-ice model. Furthermore, the presence of tensile strength in a grid cell depends on a certain coverage of a grid cell with icebergs and the present sea-ice concentration in this cell. The used threshold $A_{\text{icebergs}}<\mathit{\pi }\frac{(\mathrm{0.5}\sqrt{\left|K\right|}{)}^{\mathrm{2}}}{\left|K\right|}$ is motivated by the grounded iceberg test case and should be evaluated in more realistic setups. This also holds true for the linear dependency of the tensile strength on the sea-ice concentration. The functional relation as well as the choice of the parameter c_tensile needs to be further investigated in the context of observations.

Another assumption that could be relaxed in our ice-mélange model is the usage of uniform drag coefficients for ocean and air. Instead, different values for sea ice and icebergs could be applied. So far, the applied particle realization has not represented mechanical breakup of an iceberg or mechanical bonding of two icebergs. Implementing these processes is subject to future work.

The ice-mélange model is able to represent the formation of a polynya and the pile-up of sea ice in front of iceberg particles which are grounded on a subgrid scale. This mechanism is central for the formation of landfast sea ice in the Southern Ocean (Fraser et al., 2023). So far, the Antarctic landfast sea ice has only been poorly represented in current coupled climate models as the “fastening” mechanism due to grounded icebergs is not been taken into account. Therefore, the proposed integration of small iceberg particles into the continuum sea-ice formulation used in climate models is a promising perspective for a more realistic representation of landfast sea ice in the Southern Ocean.

6 Conclusions

We present a prototype of a dynamic hybrid ice-mélange model, which can be straightforwardly coupled to existing sea-ice components in climate models. The ice mélange is described as a joint continuum of sea ice and icebergs. While sea ice is represented as a continuum, the icebergs are modelled by particles. In order to derive a joint continuum for the ice-mélange's thickness and concentration, we integrate the particles into the sea-ice thickness and sea-ice concentration. By doing so, icebergs in the ice mélange are considered thick and compact pieces of sea ice.

Due to the use of particles in the joint continuum, we do not need to use meshes that resolve icebergs which are normally several hundred metres in size. Instead, the ice-mélange model can be simulated on the mesh resolution used for sea-ice components in climate models. This is an appealing feature with respect to the numerical efficiency.

In the context of the hybrid ice-mélange model, sea ice is modelled based on a modification of the continuum viscous–plastic sea-ice rheology, which is currently the most used material law for sea ice in climate models. Icebergs are introduced into the viscous–plastic rheology by a strength parameterization, which is used in order to prevent icebergs from diffusing.

The hybrid model is validated through a series of idealized setups that represent situations observed in nature. The setups demonstrate that the integrated icebergs, represented as thick and compact pieces of sea ice, maintain their shape under high pressure or difficult wind conditions due to the strength modification. Furthermore, we show that the hybrid ice-mélange model is capable of simulating a polynya due to subgrid iceberg grounding. These examples highlight situations where this modelling framework is beneficial. These are setups where the sea-ice cover is dense and the geometry of the fjord is complex such that sea-ice–iceberg interactions are important.

In conclusion, the use of particles and the simple extension of the viscous–plastic sea-ice material law make the hybrid model a promising approach to efficiently integrate ice mélange into climate models.