Updated October 2026: this post has been rewritten with new, more robust code.
Prompted by this tweet and its campaign for icebergs to be drawn in their most stable orientation, here is a Python script modelling the dynamics of a two-dimensional iceberg, which starts in some arbitrary position and orientation and relaxes, under gravitational and buoyant forces, to its most stable configuration.
Icebergs in adverts and cartoons are usually tall and pointed, standing upright with a small tip poking above the water and a deep keel below it. Real icebergs rarely float like that, just as a cork floats "on its side", with its longest axis parallel to the water's surface rather than bobbing around vertically.
The model
The iceberg is a rigid body, a polygon representing a cross-section through a long iceberg, so its area, \(A\), is really a volume per metre of its length, and so on. Two forces act on it:
- its weight, \(\boldsymbol{F}_\mathrm{g} = -Mg\,\hat{\boldsymbol{\jmath}}\), acting at its centre of mass, \(G\), where \(M = \rho_\mathrm{ice}A\) is its mass and \(g\) the acceleration due to gravity;
- the buoyant force, \(\boldsymbol{F}_\mathrm{b} = \rho_\mathrm{w}A_\mathrm{sub}g\,\hat{\boldsymbol{\jmath}}\), which by Archimedes' principle is equal to the weight of the water displaced by the submerged part of the iceberg (of area \(A_\mathrm{sub}\)). It acts at the centre of buoyancy, \(B\): the centroid of the submerged part.
The sea's surface is the line \(y = 0\), and the densities of ice and sea water are taken to be \(\rho_\mathrm{ice} = 0.9\;\mathrm{t\,m^{-3}}\) and \(\rho_\mathrm{w} = 1.027\;\mathrm{t\,m^{-3}}\).
The sum of the two forces accelerates the centre of mass and, where they are not in line, they exert a torque about \(G\), \(\tau = (\boldsymbol{r}_B - \boldsymbol{r}_G) \times \boldsymbol{F}_\mathrm{b}\), which turns the iceberg. Since \(\boldsymbol{F}_\mathrm{b}\) is vertical, only the horizontal distance between \(B\) and \(G\) matters, and the equations of motion are
where \(\theta\) is the angle through which the iceberg has turned (anticlockwise) and \(I\) is its moment of inertia about \(G\). The last terms are a crude model of the water's resistance to motion: linear damping at rates \(\gamma_\mathrm{v} = 5\;\mathrm{s^{-1}}\) and \(\gamma_\mathrm{r} = 1\;\mathrm{s^{-1}}\), in proportion to the fraction of the iceberg under water, \(f = A_\mathrm{sub}/A\).
The iceberg is in equilibrium when the forces balance, so that a fraction \(\rho_\mathrm{ice}/\rho_\mathrm{w} = 0.876\) of it is under water, and when \(B\) is directly below (or above) \(G\), so that there is no torque.
Stability
An equilibrium is not necessarily stable, however. Standing upright, the iceberg in the figure above has \(B\) directly below \(G\) (by 0.75 m), but when it tilts slightly, the submerged part changes shape and \(B\) moves sideways towards the side which has gone down. If \(B\) moves far enough, it ends up beyond \(G\) and the torque rights the iceberg; if not, the torque tips it further over.
For a small tilt, \(\theta\), \(B\) moves sideways by \(\overline{BM}\,\theta\), where the metacentric radius, \(\overline{BM} = w^3/(12A_\mathrm{sub})\), depends on the width, \(w\), of the iceberg at the waterline. The iceberg is stable if the metacentric height, \(\overline{GM} = \overline{BM} - \overline{BG}\), is positive. Our upright iceberg is only \(w = 6.0\;\mathrm{m}\) wide at the waterline, so \(\overline{BM} = 0.28\;\mathrm{m}\) and \(\overline{GM} = -0.47\;\mathrm{m}\): it is unstable, and topples over.
Show derivationHide derivation
Tilting the iceberg by a small angle \(\theta\) about the centre of its waterline lifts a thin triangular wedge of ice out of the water on one side and pushes a similar wedge in on the other. Each wedge has legs \(w/2\) and \((w/2)\theta\), and so area \(w^2\theta/8\), and its centroid is a distance \(w/3\) from the centre of the waterline.
The submerged area is unchanged, but an area \(w^2\theta/8\) has moved a distance \(2w/3\) across, so its centroid, \(B\), moves sideways by
The torque on the iceberg is then \((\overline{BM} - \overline{BG})\,\theta\,Mg = \overline{GM}\,\theta\,Mg\) (taking \(\overline{BG}\) as positive when \(B\) is below \(G\)), which turns it back upright if \(\overline{GM} > 0\).
The code
The polygon's area, centroid and moment of inertia are calculated from its vertex coordinates (see below). The submerged part of the iceberg is found by clipping the polygon to the region \(y \le 0\) (a single step of the Sutherland–Hodgman algorithm): walking around its edges, keeping the vertices below the waterline, and adding a vertex wherever an edge crosses it. This works for any polygon, including concave ones which cross the waterline several times.
The state of the iceberg is the vector \((x_G, y_G, \theta, \dot{x}_G, \dot{y}_G, \dot{\theta})\), and the equations of motion are integrated with SciPy's solve_ivp. The buoyant force changes abruptly as an iceberg dropped from a height hits the water, so the integrator is not allowed to take steps longer than one frame of the animation (0.1 s).
For a polygon with vertices \((x_i, y_i)\), \(i = 0, 1, \ldots, N-1\), ordered anticlockwise and with \((x_N, y_N) \equiv (x_0, y_0)\), let \(c_i = x_iy_{i+1} - x_{i+1}y_i\). Then the area (the "shoelace formula") is
its centroid is at
and its polar second moment of area about the origin is
For an iceberg of uniform density, with its centre of mass at the origin, the moment of inertia is \(I = \rho_\mathrm{ice}J\).
To run the script, which needs NumPy, SciPy and Matplotlib:
python floating_iceberg.py # an irregular iceberg dropped into the sea
python floating_iceberg.py --shape kite # the advert iceberg
python floating_iceberg.py --shape kite kite.mp4 # save the animation instead
To model your own iceberg, add an entry to the SHAPES dictionary: its vertices (in metres), its initial angle (in radians), the initial position of its centre of mass (or None to start it floating at its waterline) and the duration of the simulation. The damping rates are arguments of the Iceberg class.
"""Animate an iceberg in the sea, finding its equilibrium orientation.
Choose the iceberg's shape with --shape:
irregular: an irregular iceberg dropped into the sea from a height;
kite: the tall, pointed iceberg of adverts and cartoons, released floating
upright with only its tip showing, and nudged slightly.
Run without a filename to show the animation, or pass one (e.g. iceberg.mp4)
to save it instead.
The model is a cross-section of a long iceberg, so every area is really a
volume per unit length and every "area density" a mass per unit length: with
lengths in m and densities in t.m-3, masses are in tonnes and forces in kN,
per metre of the iceberg's length. The sea surface is the line y = 0.
The iceberg's state is the vector y = (X, Y, theta, U, V, omega): (X, Y) is
the position of its centre of mass, G; theta is its angle of rotation,
anticlockwise, from its initial shape; (U, V) is the velocity of G and
omega = dtheta/dt its angular velocity. Gravity acts at G and buoyancy at the
centre of buoyancy, B: the centroid of the submerged part of the iceberg. The
water's resistance to motion is modelled crudely as linear damping, in
proportion to how much of the iceberg is under water.
A polygon is an (N,2) array of its vertex coordinates, (x, y), without the
first vertex repeated at the end.
"""
import argparse
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from matplotlib.patches import FancyArrowPatch
from scipy.integrate import solve_ivp
from scipy.optimize import brentq
# Acceleration due to gravity, m.s-2.
g = 9.81
# Densities of sea water and ice, t.m-3.
RHO_SEA, RHO_ICE = 1.027, 0.9
SEA_COLOUR = '#8888ff'
ICE_COLOUR = '#ddddff'
SUBMERGED_ICE_COLOUR = '#ffdddd'
G_COLOUR, B_COLOUR = 'r', 'b'
# Our two-dimensional icebergs, each defined as a polygon (m) with its
# initial angle (rad) and the initial position of its centre of mass: G0 is
# None to start it floating at its equilibrium waterline, and t_max (s) is
# the duration of the simulation.
SHAPES = {
'irregular': {
'poly': [(3,0), (3,3), (1,5), (4,7), (0,12), (1,15), (4,17), (6,14),
(7,14), (8,12), (7,10), (7,7), (5,1)],
'theta0': -0.2, 'G0': (6, 15), 't_max': 25,
},
'kite': {
'poly': [(0,6), (-4,2), (0,-12), (4,2)],
'theta0': 0.5, 'G0': None, 't_max': 35,
},
}
# The time between animation frames (s).
dt = 0.1
def _cross_terms(xy):
"""Return x_i, y_i, x_{i+1}, y_{i+1} and x_i.y_{i+1} - x_{i+1}.y_i."""
x, y = xy.T
xn, yn = np.roll(x, -1), np.roll(y, -1)
return x, y, xn, yn, x * yn - xn * y
def area(xy):
"""Return the signed area of polygon xy (the "shoelace" formula).
The area is positive if the vertices are ordered anticlockwise. A
polygon with fewer than three vertices has zero area.
"""
if len(xy) < 3:
return 0.0
*_, c = _cross_terms(xy)
return np.sum(c) / 2
def centroid(xy):
"""Return the centroid (centre of mass at uniform density) of polygon xy."""
x, y, xn, yn, c = _cross_terms(xy)
return np.array((np.sum((x + xn) * c), np.sum((y + yn) * c))) / 6 / area(xy)
def second_moment(xy):
"""Return the polar second moment of area of polygon xy about the origin.
The polygon is assumed to be ordered anticlockwise.
"""
x, y, xn, yn, c = _cross_terms(xy)
Ix = np.sum(c * (y**2 + y * yn + yn**2)) / 12
Iy = np.sum(c * (x**2 + x * xn + xn**2)) / 12
# Perpendicular axis theorem.
return Ix + Iy
def transform(xy, theta, X=(0, 0)):
"""Rotate polygon xy by theta about the origin, then translate it by X."""
s, c = np.sin(theta), np.cos(theta)
return xy @ np.array(((c, s), (-s, c))) + np.asarray(X)
def clip_below(xy, level=0.0):
"""Return the part of polygon xy which lies at or below y = level.
This is the Sutherland-Hodgman algorithm for a single half-plane. If xy
is concave and crosses the line more than twice, the pieces below it are
returned as one polygon joined by zero-width "bridges" along y = level:
these do not change its area or centroid. If no part of xy is below the
line, an empty (0,2) array is returned.
"""
out = []
for P, Q in zip(xy, np.roll(xy, -1, axis=0)):
P_below, Q_below = P[1] <= level, Q[1] <= level
if P_below:
out.append(P)
if P_below != Q_below:
# This edge crosses the line: add the crossing point.
out.append(P + (level - P[1]) / (Q[1] - P[1]) * (Q - P))
return np.array(out).reshape(-1, 2)
def cross2(r, F):
"""Return the z-component of the cross product of 2D vectors r and F."""
return r[0] * F[1] - r[1] * F[0]
class Iceberg:
"""A rigid polygonal iceberg of uniform density."""
def __init__(self, poly, rho=RHO_ICE, rho_water=RHO_SEA, heave_damping=5,
roll_damping=1):
"""Initialize the iceberg.
poly is a sequence of (x, y) vertices for the iceberg's shape and rho
is its density. The damping rates (s-1) apply to the linear and
angular velocities, respectively, when it is fully submerged.
"""
xy = np.asarray(poly, dtype=float)
if area(xy) < 0:
# Order the vertices anticlockwise.
xy = xy[::-1]
self.A = area(xy)
self.M = rho * self.A
# The iceberg's shape in its own frame of reference, with its centre
# of mass at the origin.
self.xy0 = xy - centroid(xy)
self.Iz = rho * second_moment(self.xy0)
self.rho, self.rho_water = rho, rho_water
self.heave_damping, self.roll_damping = heave_damping, roll_damping
def world(self, X, theta):
"""Return the iceberg's polygon with G at X, rotated by theta."""
return transform(self.xy0, theta, X)
def hydrostatics(self, X, theta):
"""Return the submerged polygon, its area and the centre of buoyancy.
If the iceberg is out of the water, the centre of buoyancy is None.
"""
submerged = clip_below(self.world(X, theta))
A_sub = area(submerged)
B = centroid(submerged) if A_sub > 0 else None
return submerged, A_sub, B
def equilibrium_height(self, theta):
"""Return the height of G at which the iceberg floats at angle theta.
At this height the buoyant force balances the iceberg's weight, but
there is a torque on it unless B lies directly below (or above) G.
"""
A_float = self.M / self.rho_water
y = transform(self.xy0, theta)[:, 1]
def excess_area(Y):
return area(clip_below(self.world((0, Y), theta))) - A_float
return brentq(excess_area, -y.max(), -y.min())
def rhs(self, t, y):
"""Return dy/dt for the state vector y, for scipy's solve_ivp."""
X, Y, theta, U, V, omega = y
G = np.array((X, Y))
_, A_sub, B = self.hydrostatics(G, theta)
Fg = np.array((0, -self.M * g))
Fb = np.array((0, self.rho_water * A_sub * g))
tau = cross2(B - G, Fb) if B is not None else 0
a = (Fg + Fb) / self.M
alpha = tau / self.Iz
# Damp the motion in proportion to the submerged fraction.
f_sub = A_sub / self.A
a -= self.heave_damping * f_sub * np.array((U, V))
alpha -= self.roll_damping * f_sub * omega
return np.array((U, V, omega, a[0], a[1], alpha))
def simulate(iceberg, G0, theta0, t_max):
"""Integrate the motion of iceberg from rest with G at G0, angle theta0."""
y0 = (*G0, theta0, 0, 0, 0)
t = np.arange(0, t_max + dt / 2, dt)
# The buoyant force changes abruptly as the iceberg enters the water, so
# don't let the integrator take steps longer than a frame.
sol = solve_ivp(iceberg.rhs, (0, t_max), y0, t_eval=t, max_step=dt)
return sol.t, sol.y.T
def animate(iceberg, t, states):
"""Return a FuncAnimation of the iceberg following the states at times t."""
# Fix the plot limits so that the whole trajectory fits.
all_xy = np.vstack([iceberg.world(s[:2], s[2]) for s in states])
(xmin, ymin), (xmax, ymax) = all_xy.min(axis=0) - 2, all_xy.max(axis=0) + 2
# Size the figure to fit the plot: 6 inches tall, with room for labels.
width = 6 * (xmax - xmin) / (ymax - ymin) + 1
fig, ax = plt.subplots(figsize=(width, 6.5), layout='constrained')
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
# Equal distances in the x- and y-directions should look the same.
ax.set_aspect('equal')
ax.set_xlabel('x / m')
ax.set_ylabel('y / m')
# The sea! The sea!
ax.add_patch(plt.Rectangle((xmin, ymin), xmax - xmin, -ymin,
fc=SEA_COLOUR))
ice_patch = ax.add_patch(plt.Polygon(iceberg.xy0, fc=ICE_COLOUR, ec='k'))
submerged_patch = ax.add_patch(plt.Polygon(iceberg.xy0,
fc=SUBMERGED_ICE_COLOUR, ec='k'))
# The centre of mass, G (red), with the force of gravity, F_g, acting
# on it, and the centre of buoyancy, B (blue), with the buoyant force,
# F_b, acting on it. Each arrow's length is proportional to its force:
# the weight's arrow is half as long as the square root of the
# iceberg's area.
force_scale = 0.5 * np.sqrt(iceberg.A) / (iceberg.M * g)
def arrow(colour):
return ax.add_patch(FancyArrowPatch((0, 0), (0, 1), color=colour,
arrowstyle='-|>', mutation_scale=15, lw=1.5,
shrinkA=0, shrinkB=0, zorder=3))
def label(text, colour, ha):
return ax.text(0, 0, text, color=colour, ha=ha, va='center',
fontsize=14, zorder=3)
Fg_arrow, Fb_arrow = arrow(G_COLOUR), arrow(B_COLOUR)
Fg_label = label(r'$\mathbfit{F}_\mathrm{g}$', G_COLOUR, 'center')
Fb_label = label(r'$\mathbfit{F}_\mathrm{b}$', B_COLOUR, 'center')
G_patch = ax.add_patch(plt.Circle((0, 0), 0.2, fc=G_COLOUR, zorder=4))
B_patch = ax.add_patch(plt.Circle((0, 0), 0.2, fc=B_COLOUR, zorder=4))
G_label = label(r'$G$', G_COLOUR, 'left')
B_label = label(r'$B$', B_COLOUR, 'right')
time_text = ax.text(0.02, 0.95, '', transform=ax.transAxes)
def update(i):
G, theta = states[i, :2], states[i, 2]
ice_patch.set_xy(iceberg.world(G, theta))
submerged, A_sub, B = iceberg.hydrostatics(G, theta)
in_water = B is not None
for artist in (submerged_patch, B_patch, B_label, Fb_arrow, Fb_label):
artist.set_visible(in_water)
Fg_end = G - (0, force_scale * iceberg.M * g)
Fg_arrow.set_positions(G, Fg_end)
Fg_label.set_position(Fg_end - (0, 1))
G_patch.center = G
G_label.set_position(G + (0.4, 0))
if in_water:
submerged_patch.set_xy(submerged)
Fb_end = B + (0, force_scale * iceberg.rho_water * A_sub * g)
Fb_arrow.set_positions(B, Fb_end)
Fb_label.set_position(Fb_end + (0, 1))
B_patch.center = B
B_label.set_position(B - (0.4, 0))
time_text.set_text(f't = {t[i]:.1f} s')
return (ice_patch, submerged_patch, Fg_arrow, Fb_arrow, Fg_label,
Fb_label, G_patch, B_patch, G_label, B_label, time_text)
return animation.FuncAnimation(fig, update, len(t), blit=True,
interval=dt * 1000)
if __name__ == '__main__':
parser = argparse.ArgumentParser(
description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument('filename', nargs='?',
help='save the animation to this file')
parser.add_argument('--shape', choices=SHAPES, default='irregular',
help='the shape of the iceberg')
args = parser.parse_args()
shape = SHAPES[args.shape]
iceberg = Iceberg(shape['poly'])
theta0, G0 = shape['theta0'], shape['G0']
if G0 is None:
G0 = (0, iceberg.equilibrium_height(theta0))
t, states = simulate(iceberg, G0, theta0, shape['t_max'])
ani = animate(iceberg, t, states)
if args.filename:
ani.save(args.filename)
else:
plt.show()
Limitations
This is a two-dimensional model of an iceberg whose cross-section is the same along its length, and the water is treated very simply: there are no waves, and the water's resistance to motion is modelled only crudely by damping. Real icebergs also melt and break up, which changes their shape and can suddenly make them unstable: an iceberg which capsizes as a result can release a great deal of energy. See, for example, the article on the use of catastrophe theory to analyse the stability and toppling of icebergs in the Annals of Glaciology.