Module 8: Climate Change & Conservation — A Species at the Edge

Every element of emperor penguin biology that the preceding modules have examined — the 120-day fast, the two-voice syrinx, the 30-km/h hydrodynamic burst, the huddle thermodynamics — is tuned to a specific climatic regime: the annual fast-ice cycle of the Southern Ocean. That regime is changing. Stephanie Jenouvrier and colleagues (2014, Nature Climate Change; 2021, Global Change Biology) combined CMIP5/CMIP6 sea-ice projections with matrix-population demography to predict a 70% global decline of emperor populations by 2100 under RCP 8.5, with 50 of 54 known colonies at high risk of quasi-extinction. Fretwell & Trathan’s (2019) Halley Bay collapse, the ongoing H5N1 panzootic (Banyard 2022; Ashcroft 2024), krill decline (Atkinson 2019), and increasing ecotourism together compound the pressure. This module reviews the demographic modelling, the Antarctic governance framework, and the conservation measures that may slow but not reverse the trajectory.

1. The Jenouvrier Projections (2014, 2021)

Jenouvrier, Caswell, Barbraud, Holland, Strøeve & Weimerskirch (2014, Nature Climate Change) published the first rigorous climate-coupled demographic projection for emperor penguins. Their framework was a coupled Leslie matrix population model driven by CMIP5 sea-ice extent (SIE) projections under RCP 2.6, 4.5, and 8.5, with vital rates modulated by SIE anomalies estimated from the 50-year Pointe Géologie mark-recapture dataset.

\[ \mathbf{N}(t+1) \;=\; \mathbf{L}[\text{SIE}(t)] \cdot \mathbf{N}(t) \]

Discrete-time projection: \(\mathbf{N}\) is the age-structured population vector; \(\mathbf{L}\) is the climate-dependent Leslie matrix with SIE-modulated fecundity and survivorship.

Key findings (Jenouvrier 2014, 2021):

Under RCP 8.5 / SSP 5-8.5: global population decline of 81% (CI 55–95%) by 2100.
Under RCP 4.5 / SSP 2-4.5: decline of 31% (CI 15–54%).
Under RCP 2.6 / SSP 1-2.6: decline of 19% (CI −10–50%) — the only scenario compatible with species persistence.
50 of 54 known colonies project quasi-extinction by 2100 under RCP 8.5.
Colonies at higher latitudes (deeper ice) slightly buffer the decline; low-latitude colonies (e.g., Pointe Géologie, Halley Bay) are most exposed.

The 2021 update incorporated Antarctic sea-ice dynamics rather than only extent: earlier break-up, later formation, and more variable fast-ice thickness. The updated projection showed additional risk from intra-seasonal ice variability — even when annual mean SIE is preserved, a break-up event during chick-rearing can cause 100% breeding failure at an affected colony (as observed at Halley Bay).

2. Life-History Coupling to Fast Ice

The emperor breeding cycle is locked to the annual fast-ice schedule in a way not replicated by any other seabird. The critical couplings (Barbraud & Weimerskirch 2001; Trathan et al. 2011):

Colony formation (March–April): requires fast-ice formation within 10 km of a foraging ocean and stable through June.
Incubation (May–July): requires continuous fast-ice support of the colony for ~70 days. Ice break-up during incubation = 100% egg loss.
Chick-rearing (August–November): requires fast-ice stability until chicks can swim (~day 150).
Fledging (December): timed to coincide with ice break-up opening swim routes to prey.

Too little ice kills the colony; too much ice (a historical rarity, now exceptional) imposes energetic cost by extending the walking distance from ice edge to foraging area. The optimal ice regime — the Goldilocks zone — is narrow, and climate projections show that zone shifting poleward at a rate faster than colonies can relocate.

\[ \phi(\text{SIE}) \;=\; \exp\!\left(-\frac{(\text{SIE} - \text{SIE}^{*})^2}{2\sigma_{\text{SIE}}^2}\right) \]

Vital-rate multiplier: Gaussian bell centred on the historical optimum\(\text{SIE}^{*}\), with width \(\sigma_{\text{SIE}}\).

3. The Halley Bay Collapse (2016–2019)

Fretwell & Trathan (2019, Antarctic Science) documented the near-total collapse of the Halley Bay colony in the Weddell Sea — one of the five largest emperor colonies globally. Using Sentinel-2 and Landsat satellite imagery, they showed that fast-ice break-up during the 2016 austral spring caused complete breeding failure (estimated 14 000–25 000 chick deaths). The colony did not recover in subsequent years; much of the surviving breeding population relocated to the nearby Dawson-Lambton colony.

The Halley Bay event was the first directly-observed climate-driven colony collapse at this scale. Its causes are instructive:

An abnormally strong early-austral-spring windstorm (October 2016) broke up fast-ice adjacent to the colony.
Chicks lacked swim-ready feathers (moult completes ~December); the sudden exposure to open water was lethal.
Three successive years of fast-ice failure followed (2017–2019), preventing recolonisation of the original site.
Breeding adults dispersed to alternative colonies. Whether these “climate refugees” successfully integrated into those populations remains uncertain.

The event showed that emperor colonies are vulnerable not only to gradual SIE trends but to stochastic ice-break-up events. Even if mean SIE is preserved, increased interannual variability can cascade through the colony system via repeated breeding failures.

4. Schematic: Compressed Breeding Window

Simulation 1: CMIP6 Climate Forcing & Leslie-Matrix Population Model

Build a 6-age-class Leslie matrix with SIE-dependent vital rates. Generate ensembles of stochastic SIE trajectories 2020–2100 for SSP 1-2.6, SSP 2-4.5, and SSP 5-8.5. Run 200-member Monte Carlo demographic simulations. Quantify the probability of quasi-extinction (adults < 10% of 2020 level) by 2100.

Python

script.py224 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Climate-driven demographic projections of emperor penguin colonies.
#
# Based on Jenouvrier et al. (2014, Nature Climate Change; 2021, Global Change
# Biology) who coupled CMIP5/CMIP6 sea-ice projections to a Leslie-matrix
# population model for emperor penguins. We construct a compact version of the
# method:
#
#   1. Prescribe SSP 1-2.6, SSP 2-4.5, SSP 5-8.5 sea-ice trajectories for
#      Antarctic sea-ice extent (SIE) in the Pointe Geologie sector, scaled to
#      match Jenouvrier's published decadal anomalies.
#   2. Map SIE to a vital-rate multiplier phi(SIE) relating fast-ice stability
#      to breeding success (Barbraud & Weimerskirch 2001).
#   3. Build a 6-age-class Leslie matrix L(t) with SIE-dependent fecundity and
#      survivorship.
#   4. Iterate 2020-2100 on an ensemble of stochastic realisations with
#      interannual SIE noise to get colony-level extinction probability.

np.random.seed(42)

# Time horizon
years = np.arange(2020, 2101)
n_y = len(years)

# Sea-ice baseline and CMIP-scenario sensitivities (scaled from Jenouvrier 2021)
SIE_2020 = 1.0       # normalised 2020 baseline (Pointe Geologie sector)

# Linear+noise trends by scenario (% change per decade)
trends = {
    'SSP 1-2.6': -0.015,   # -1.5%/decade
    'SSP 2-4.5': -0.040,   # -4.0%/decade
    'SSP 5-8.5': -0.080,   # -8.0%/decade
}
noise_sigma = 0.06  # year-to-year noise

def generate_SIE(trend_dec, seed):
    rng = np.random.default_rng(seed)
    years_since_2020 = years - 2020
    sie = SIE_2020 * (1.0 + trend_dec * (years_since_2020 / 10.0))
    sie = sie + noise_sigma * rng.standard_normal(n_y)
    return np.clip(sie, 0.05, 1.4)

# Vital-rate mapping: phi(SIE) in [0, 1], Gaussian-bell shape
# Too little fast-ice -> collapse; too much -> energy cost of longer ice walk
# Peak at SIE = 1.05 (historical optimum), width 0.25
def phi_SIE(sie):
    return np.exp(-((sie - 1.05) / 0.25) ** 2)

# Age classes 0,1,2,3,4,5+; age 5+ is the breeding class
# Baseline Leslie matrix
# Fecundity f = 0.5 * phi * breeding_success (baseline 0.6)
# Adult survival s_adult = 0.95 * (0.6 + 0.4*phi)
# Juvenile survival s_juv = 0.10 * (0.5 + 0.5*phi)
def build_L(sie):
    p = phi_SIE(sie)
    s_adult = 0.95 * (0.6 + 0.4 * p)
    s_juv = 0.1 * (0.5 + 0.5 * p)
    fec = 0.25 * p * 0.6  # per-capita 5+ fecundity (females only)
    L = np.zeros((6, 6))
    L[0, 5] = fec
    L[1, 0] = s_juv
    L[2, 1] = 0.82
    L[3, 2] = 0.85
    L[4, 3] = 0.90
    L[5, 4] = 0.92
    L[5, 5] = s_adult
    return L

# Initial population distribution (Pointe Geologie, nominal 10 000 pairs =
# 20 000 adults)
N0 = np.array([2500, 2000, 1800, 1600, 1500, 10600], dtype=float)

def simulate_scenario(trend_dec, seed):
    sie = generate_SIE(trend_dec, seed)
    N = N0.copy()
    totals = np.zeros(n_y)
    totals[0] = N.sum()
    adults = np.zeros(n_y)
    adults[0] = N[5]
    sie_out = sie.copy()
    for i in range(1, n_y):
        L = build_L(sie[i])
        N = L @ N
        # demographic stochasticity
        N = N + 0.02 * np.sqrt(np.maximum(N, 0)) * np.random.randn(6)
        N = np.maximum(N, 0.0)
        totals[i] = N.sum()
        adults[i] = N[5]
    return totals, adults, sie_out

# Monte Carlo ensemble: 200 realisations per scenario, track 2100 population
N_ENS = 200
ensembles = {}
extinction = {}
for name, trend in trends.items():
    totals_all = np.zeros((N_ENS, n_y))
    adults_all = np.zeros((N_ENS, n_y))
    sie_all = np.zeros((N_ENS, n_y))
    for k in range(N_ENS):
        totals, adults, sie = simulate_scenario(trend, seed=k + hash(name) % 9999)
        totals_all[k] = totals
        adults_all[k] = adults
        sie_all[k] = sie
    ensembles[name] = (totals_all, adults_all, sie_all)
    # Quasi-extinction: adult count below 10% of N0[5]
    adult_2100 = adults_all[:, -1]
    ext_prob = float(np.mean(adult_2100 < 0.10 * N0[5]))
    extinction[name] = ext_prob

# Figure
fig, axes = plt.subplots(2, 3, figsize=(17, 10))
fig.patch.set_facecolor('#0a0a1a')
for ax in axes.flat:
    ax.set_facecolor('#0a0a1a')
    ax.tick_params(colors='#e2e8f0')
    for spine in ax.spines.values():
        spine.set_color('#475569')
    ax.grid(True, color='#334155', alpha=0.35)

colours = {'SSP 1-2.6': '#a3e635', 'SSP 2-4.5': '#fbbf24', 'SSP 5-8.5': '#f43f5e'}

ax = axes[0, 0]
for name, (_, _, sie_all) in ensembles.items():
    med = np.median(sie_all, axis=0)
    q05, q95 = np.quantile(sie_all, [0.05, 0.95], axis=0)
    ax.fill_between(years, q05, q95, color=colours[name], alpha=0.2)
    ax.plot(years, med, color=colours[name], lw=2.2, label=name)
ax.set_xlabel('Year', color='#e2e8f0')
ax.set_ylabel('Normalised SIE (fast ice)', color='#e2e8f0')
ax.set_title('Sea-ice extent projection (CMIP6 scenarios)',
             color='#fde68a', fontsize=12, fontweight='bold')
ax.legend(facecolor='#1e293b', edgecolor='#475569', labelcolor='#e2e8f0')

ax = axes[0, 1]
sie_grid = np.linspace(0.3, 1.4, 200)
ax.plot(sie_grid, phi_SIE(sie_grid), color='#60a5fa', lw=2.4)
ax.axvline(1.0, color='#fde68a', ls=':', alpha=0.7, label='2020 baseline')
ax.axvline(1.05, color='#a3e635', ls='--', alpha=0.7, label='optimum')
ax.set_xlabel('Normalised SIE', color='#e2e8f0')
ax.set_ylabel('Vital-rate multiplier phi', color='#e2e8f0')
ax.set_title('Fast-ice -> vital-rate mapping',
             color='#fde68a', fontsize=12, fontweight='bold')
ax.legend(facecolor='#1e293b', edgecolor='#475569', labelcolor='#e2e8f0')

ax = axes[0, 2]
for name, (_, adults_all, _) in ensembles.items():
    med = np.median(adults_all, axis=0)
    q05, q95 = np.quantile(adults_all, [0.05, 0.95], axis=0)
    ax.fill_between(years, q05, q95, color=colours[name], alpha=0.2)
    ax.plot(years, med, color=colours[name], lw=2.2, label=name)
ax.axhline(0.10 * N0[5], color='#94a3b8', ls=':', alpha=0.6,
           label='quasi-extinction')
ax.set_xlabel('Year', color='#e2e8f0')
ax.set_ylabel('Breeding adults', color='#e2e8f0')
ax.set_title('Adult population trajectory (10 000-pair colony)',
             color='#fde68a', fontsize=12, fontweight='bold')
ax.legend(facecolor='#1e293b', edgecolor='#475569', labelcolor='#e2e8f0', fontsize=9)

ax = axes[1, 0]
for name, (totals_all, _, _) in ensembles.items():
    # Express as % of initial N
    pct = 100.0 * totals_all / N0.sum()
    med = np.median(pct, axis=0)
    q05, q95 = np.quantile(pct, [0.05, 0.95], axis=0)
    ax.fill_between(years, q05, q95, color=colours[name], alpha=0.2)
    ax.plot(years, med, color=colours[name], lw=2.2, label=name)
ax.axhline(30.0, color='#f43f5e', ls='--', alpha=0.6, label='70% decline')
ax.set_xlabel('Year', color='#e2e8f0')
ax.set_ylabel('Colony (% of 2020 baseline)', color='#e2e8f0')
ax.set_title('Jenouvrier 2014: 70% decline by 2100 under SSP 5-8.5',
             color='#fde68a', fontsize=12, fontweight='bold')
ax.legend(facecolor='#1e293b', edgecolor='#475569', labelcolor='#e2e8f0', fontsize=9)

ax = axes[1, 1]
names = list(extinction.keys())
vals = [extinction[n] * 100.0 for n in names]
bar_colors = [colours[n] for n in names]
ax.bar(names, vals, color=bar_colors, edgecolor='#fde68a', lw=1.5)
for i, v in enumerate(vals):
    ax.text(i, v + 2.0, f'{v:.0f}%', ha='center', color='#fde68a', fontsize=12)
ax.set_ylabel('Quasi-extinction probability (%)', color='#e2e8f0')
ax.set_title('Quasi-extinction by 2100 (200-member ensemble)',
             color='#fde68a', fontsize=12, fontweight='bold')
ax.set_ylim(0, 105)

ax = axes[1, 2]
# Sensitivity: how does trend slope map to 2100 population?
trend_sweep = np.linspace(-0.12, 0.0, 25)
pop_2100 = []
for tr in trend_sweep:
    sample = []
    for k in range(30):
        totals, _, _ = simulate_scenario(tr, seed=k)
        sample.append(totals[-1])
    pop_2100.append(np.mean(sample))
ax.plot(np.array(trend_sweep) * 100.0, np.array(pop_2100) / N0.sum() * 100.0,
        color='#60a5fa', lw=2.4)
ax.scatter(np.array([t for t in trends.values()]) * 100.0,
           [100.0 * ensembles[n][0][:, -1].mean() / N0.sum() for n in trends],
           c=[colours[n] for n in trends], s=80, edgecolor='#fde68a', zorder=5)
ax.axhline(30.0, color='#f43f5e', ls='--', alpha=0.6)
ax.set_xlabel('SIE decadal trend (%)', color='#e2e8f0')
ax.set_ylabel('2100 population (% of 2020)', color='#e2e8f0')
ax.set_title('Sensitivity to SIE trend',
             color='#fde68a', fontsize=12, fontweight='bold')

plt.tight_layout()
plt.savefig('output.png', dpi=120, bbox_inches='tight', facecolor='#0a0a1a')

for name in trends:
    pct = 100.0 * ensembles[name][0][:, -1].mean() / N0.sum()
    adult_final = ensembles[name][1][:, -1].mean()
    print('%s: 2100 colony mean = %.1f%% of 2020, mean adults = %.0f, P(quasi-ext) = %.0f%%'
          % (name, pct, adult_final, extinction[name] * 100.0))
print('Integrated population-years lost 2020-2100 SSP 5-8.5 vs SSP 1-2.6 (%):')
pop_hi = ensembles['SSP 5-8.5'][0].mean(axis=0)
pop_lo = ensembles['SSP 1-2.6'][0].mean(axis=0)
integrated = float(np.trapezoid(pop_lo - pop_hi, years))
print('  integrated = %.2e pop-years' % integrated)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

5. Sectoral Variation: Ross Sea, Weddell Sea, East Antarctica

Antarctic sea ice does not behave as a uniform field. Regional circulation patterns and atmospheric teleconnections impose sharp sector-level heterogeneity (Parkinson 2019, PNAS):

Ross Sea: historical SIE increase 1978–2015, but reversed from 2016 onward. Colonies here (e.g., Cape Washington, Cape Roget) retained stable populations until recently.
Weddell Sea: large interannual variability; Halley Bay collapse (2016–2019) marks the first major failure event.
Amundsen-Bellingshausen: most rapid decline in SIE; largely lacks resident emperor colonies (unsuitable bathymetry + water temperature).
East Antarctica (Adélie Land, Wilkes Land): moderate decline; Pointe Géologie colony demographics now tracked 60+ years.
Indian Ocean sector: warmer air-temperature trends; Terra Nova Bay, Mertz Glacier colonies increasingly exposed.

The sector heterogeneity implies that global projections must aggregate colony-specific vulnerabilities rather than apply a single climate-envelope. Jenouvrier (2021) built an individual-colony forcing grid at 25 km resolution, showing that 50 of 54 colonies project quasi-extinction under SSP 5-8.5 but with substantial variance in timing.

6. H5N1 Highly-Pathogenic Avian Influenza (2022–2024)

HPAI H5N1 clade 2.3.4.4b began sweeping through European seabirds in 2021 (Banyard et al. 2022). By October 2022 it had reached South Georgia; by February 2023 it was detected in South American penguin colonies; by spring 2024 confirmed cases appeared in the maritime Antarctic (Ashcroft et al. 2024).

Penguin colonies present a particular vulnerability:

Huddling: mass-action contact in dense huddles accelerates transmission. R0 estimates for H5N1 in dense penguin aggregations exceed 2.0 (Ashcroft 2024).
Naive immunity: Antarctic seabird populations lack prior exposure to clade 2.3.4.4b strains.
Case fatality: estimated 60–80% for adult penguins infected with H5N1 (based on gentoo and Adelie mortality observations).
No effective intervention: vaccination or culling strategies that work for poultry are not scalable to wild Antarctic colonies.

\[ \frac{dS}{dt} = -\beta S I / N,\quad \frac{dI}{dt} = \beta S I / N - \gamma I,\quad R_0 = \beta / \gamma \]

Classical SIR dynamics. R0 = basic reproduction number; above 1, epidemic grows. For H5N1 in penguin huddles, R0 estimated 1.5–3 depending on colony density.

Emperor-specific exposure risk through 2024 is so far lower than for Adelie or gentoo (emperors breed on sea ice rather than guano-strewn rock, and aggregate less with other seabird species). However, Weimerskirch et al. (in preparation 2025) report first confirmed H5N1 in emperor chicks at an Adélie Land colony in December 2024. Stochastic models indicate that a single introduction can cause > 50% mortality within 30 days.

Simulation 2: Stochastic SIR with Super-Spreader Events

Simulate H5N1 propagation in a 10 000-adult colony using a stochastic SEIR-D (susceptible, exposed, infectious, recovered, dead) model. Include super-spreader events with probability 2% per day generating 10–30 extra secondary cases. Run 200 realisations; compute mortality, peak-infectious, and peak-timing distributions; scan R0 sensitivity.

Python

script.py183 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Stochastic SIR model of H5N1 highly-pathogenic avian influenza propagation
# through an emperor penguin colony.
#
# Based on real outbreaks reported since 2022 in subantarctic and maritime
# Antarctic seabird colonies, and the modelling framework of Ashcroft et al.
# (2024) and Banyard et al. (2022) for HPAI-H5N1 clade 2.3.4.4b spread in
# wild-bird populations.
#
# State variables in each colony:
#   S(t) susceptible; E(t) exposed (latent); I(t) infectious; R(t) recovered;
#   D(t) dead.
# Transitions (per adult, per day):
#   S -> E at rate beta * I / N (mass-action contact in huddle)
#   E -> I at rate sigma (1/latent)
#   I -> R at rate gamma * (1 - CFR)
#   I -> D at rate gamma * CFR
#
# Super-spreader events model occasional burst contacts in dense huddle
# clusters.

from scipy.stats import gamma as scipy_gamma

np.random.seed(11)

# Parameters (Ashcroft 2024, Banyard 2022, avian influenza literature)
N0 = 10000          # initial susceptible adults
beta_base = 0.35    # contact-transmission rate per day
sigma = 0.5         # 1 / 2-day latent
gamma_rate = 0.2    # 1 / 5-day infectious
CFR = 0.60          # case fatality rate (H5N1 clade 2.3.4.4b in penguins)
days = np.arange(0, 181)

# Super-spreader mechanism: each day, small probability that the bird with
# highest contact count generates an extra 10-30 secondary cases.
def run_epidemic(seed, N0=N0, beta_base=beta_base, p_super=0.02, super_range=(10, 30)):
    rng = np.random.default_rng(seed)
    S = float(N0)
    E = 0.0
    I = 5.0  # seeded with 5 infectious adults
    R = 0.0
    D = 0.0
    traj = {'S': [S], 'E': [E], 'I': [I], 'R': [R], 'D': [D]}
    for t in days[1:]:
        N = S + E + I + R  # living birds
        # Daily contact rate modulated by huddle seasonality
        huddle_factor = 1.0 + 0.5 * np.exp(-((t - 60) / 30.0) ** 2)
        beta = beta_base * huddle_factor
        lam = beta * I / max(N, 1.0)
        new_E = rng.binomial(int(round(S)), min(1.0, lam))
        new_I = rng.binomial(int(round(E)), 1.0 - np.exp(-sigma))
        new_R = rng.binomial(int(round(I)), (1.0 - np.exp(-gamma_rate)) * (1.0 - CFR))
        new_D = rng.binomial(int(round(I)), (1.0 - np.exp(-gamma_rate)) * CFR)
        # Super-spreader event
        if rng.random() < p_super and I > 0:
            extra = rng.integers(super_range[0], super_range[1] + 1)
            extra_cases = min(extra, int(S))
            new_E += extra_cases
        S = max(0.0, S - new_E)
        E = max(0.0, E + new_E - new_I)
        I = max(0.0, I + new_I - new_R - new_D)
        R = R + new_R
        D = D + new_D
        traj['S'].append(S)
        traj['E'].append(E)
        traj['I'].append(I)
        traj['R'].append(R)
        traj['D'].append(D)
    return {k: np.array(v) for k, v in traj.items()}

# Ensemble of 200 stochastic realisations
N_ENS = 200
runs = [run_epidemic(k) for k in range(N_ENS)]

# Final mortality statistics
final_D = np.array([r['D'][-1] for r in runs])
peak_I = np.array([r['I'].max() for r in runs])
peak_day = np.array([int(np.argmax(r['I'])) for r in runs])

# R0 estimate from beta, gamma (effective)
R0 = beta_base / gamma_rate

ax = axes[0, 0]
for r in runs[:40]:
    ax.plot(days, r['I'], color='#f43f5e', alpha=0.12, lw=0.9)
mean_I = np.mean([r['I'] for r in runs], axis=0)
ax.plot(days, mean_I, color='#fbbf24', lw=2.4, label='ensemble mean')
ax.set_xlabel('Days since seed', color='#e2e8f0')
ax.set_ylabel('Infectious adults I(t)', color='#e2e8f0')
ax.set_title('Stochastic infectious trajectory (200 ensemble members)',
             color='#fde68a', fontsize=12, fontweight='bold')
ax.legend(facecolor='#1e293b', edgecolor='#475569', labelcolor='#e2e8f0')

ax = axes[0, 1]
mean_S = np.mean([r['S'] for r in runs], axis=0)
mean_R = np.mean([r['R'] for r in runs], axis=0)
mean_D = np.mean([r['D'] for r in runs], axis=0)
ax.plot(days, mean_S, color='#a3e635', lw=2.4, label='S')
ax.plot(days, mean_I, color='#f43f5e', lw=2.4, label='I')
ax.plot(days, mean_R, color='#60a5fa', lw=2.4, label='R')
ax.plot(days, mean_D, color='#fbbf24', lw=2.4, label='D')
ax.set_xlabel('Days since seed', color='#e2e8f0')
ax.set_ylabel('Population compartments', color='#e2e8f0')
ax.set_title('Mean SEIR-D compartmental dynamics',
             color='#fde68a', fontsize=12, fontweight='bold')
ax.legend(facecolor='#1e293b', edgecolor='#475569', labelcolor='#e2e8f0')

ax = axes[0, 2]
ax.hist(final_D, bins=30, color='#f43f5e', edgecolor='#fde68a', alpha=0.85)
ax.axvline(np.median(final_D), color='#fbbf24', ls='--', lw=2,
           label=f'median {np.median(final_D):.0f}')
ax.set_xlabel('Final cumulative mortality D', color='#e2e8f0')
ax.set_ylabel('Ensemble count', color='#e2e8f0')
ax.set_title('Mortality distribution (180 days)',
             color='#fde68a', fontsize=12, fontweight='bold')
ax.legend(facecolor='#1e293b', edgecolor='#475569', labelcolor='#e2e8f0')

ax = axes[1, 0]
ax.hist(peak_I, bins=30, color='#fbbf24', edgecolor='#fde68a', alpha=0.85)
ax.axvline(np.median(peak_I), color='#a3e635', ls='--', lw=2,
           label=f'median {np.median(peak_I):.0f}')
ax.set_xlabel('Peak infectious I_peak', color='#e2e8f0')
ax.set_ylabel('Ensemble count', color='#e2e8f0')
ax.set_title('Peak-infectious distribution',
             color='#fde68a', fontsize=12, fontweight='bold')
ax.legend(facecolor='#1e293b', edgecolor='#475569', labelcolor='#e2e8f0')

ax = axes[1, 1]
ax.hist(peak_day, bins=30, color='#60a5fa', edgecolor='#fde68a', alpha=0.85)
ax.axvline(np.median(peak_day), color='#fbbf24', ls='--', lw=2,
           label=f'median day {int(np.median(peak_day))}')
ax.set_xlabel('Day of peak infection', color='#e2e8f0')
ax.set_ylabel('Ensemble count', color='#e2e8f0')
ax.set_title('Peak-timing distribution',
             color='#fde68a', fontsize=12, fontweight='bold')
ax.legend(facecolor='#1e293b', edgecolor='#475569', labelcolor='#e2e8f0')

ax = axes[1, 2]
# Scan R0 by varying beta
beta_grid = np.linspace(0.10, 0.60, 12)
mortality_grid = []
for b in beta_grid:
    runs_b = [run_epidemic(k, beta_base=b) for k in range(30)]
    mortality_grid.append(np.mean([r['D'][-1] for r in runs_b]))
R0_grid = beta_grid / gamma_rate
ax.plot(R0_grid, np.array(mortality_grid) / N0 * 100.0,
        color='#a3e635', lw=2.4, marker='o', markersize=6)
ax.axvline(1.0, color='#fde68a', ls=':', alpha=0.6, label='R0 = 1')
ax.axvline(R0, color='#f43f5e', ls='--', alpha=0.6,
           label=f'baseline R0 = {R0:.1f}')
ax.set_xlabel('Basic reproduction number R0', color='#e2e8f0')
ax.set_ylabel('Final mortality (%)', color='#e2e8f0')
ax.set_title('Final mortality vs. R0',
             color='#fde68a', fontsize=12, fontweight='bold')
ax.legend(facecolor='#1e293b', edgecolor='#475569', labelcolor='#e2e8f0')

plt.tight_layout()
plt.savefig('output.png', dpi=120, bbox_inches='tight', facecolor='#0a0a1a')

print('Colony size: %d adults' % N0)
print('R0 (baseline): %.2f' % R0)
print('CFR: %.2f' % CFR)
print('Mean final mortality: %.0f (%.1f%%)' % (np.mean(final_D),
                                                 100.0 * np.mean(final_D) / N0))
print('Median peak infectious: %.0f (day %d)' % (np.median(peak_I),
                                                  int(np.median(peak_day))))
print('Integrated person-days infectious: %.2e'
      % float(np.trapezoid(mean_I, days)))

Click Run to execute the Python code

Code will be executed with Python 3 on the server

7. Krill Decline and Prey-Base Erosion

Emperors feed primarily on Antarctic silverfish (Pleuragramma antarctica) and krill (Euphausia superba) during the foraging phase. Krill abundance in the Scotia Sea has declined by 70–80% since the 1970s (Atkinson, Hill, Pakhomov, et al. 2019, Nature Climate Change). The southern shift of the krill population towards higher latitudes decouples historical prey-distribution overlap with emperor colonies.

The krill-silverfish coupling is critical during the chick-provisioning period (December–February), when adults must deliver meals every 3–10 days. Increased foraging-trip distance reduces meal frequency, increasing chick starvation risk. Le Guen et al. (2018) reported a 30% drop in meal frequency at Dumont d’Urville during El Niño years, correlating with reduced silverfish catch rates near the colony.

ENSO and SAM teleconnections

The Southern Annular Mode (SAM) and El Niño-Southern Oscillation (ENSO) drive multi-decadal variability in Antarctic sea-ice and prey distribution:

Positive SAM: strengthened westerlies, warming of the Antarctic Peninsula, southward shift of sea-ice edge, reduced krill recruitment in the Scotia Sea.
El Niño: reduces Antarctic sea-ice in the Amundsen-Bellingshausen sector; varied effect in East Antarctica.
Combined +SAM, +ENSO: compound stress on breeding success at low-latitude colonies.

8. Human Pressures: Fishing, Tourism, Pollution

Direct anthropogenic pressures on emperors have historically been modest compared to climate forcing, but they are increasing:

Krill fishery

CCAMLR (Commission for the Conservation of Antarctic Marine Living Resources) regulates krill extraction in the Southern Ocean. The precautionary catch limit for Area 48 (Atlantic sector) has stood at 620 000 tonnes/year, but concentrated effort in the Scotia Sea overlaps with key predator foraging zones. Emperor populations are not directly targeted but compete for krill-dependent silverfish.

Ecotourism

Antarctic tourism has grown from ~5 000 visitors in the 1990s to over 100 000 per austral summer by 2023. Emperor colony landings are rare but visitation to ship-accessible colonies (e.g., Snow Hill Island) can reach 2 000 visitors per season. Disturbance effects include elevated heart rate, altered breeding behaviour, and the potential for pathogen introduction from bird-to-bird via human vectors.

Pollution

Persistent organic pollutants (POPs), mercury, and microplastics accumulate through the Antarctic food web. Bustamante et al. (2021) reported Hg concentrations 20% higher in modern emperor feathers compared to 19th-century museum specimens. Microplastic contamination has been detected in penguin guano at multiple colonies.

9. Remote Sensing of Colonies (Sentinel-2, Landsat, Maxar)

Emperor colonies are logistically inaccessible through most of the year. Remote sensing has revolutionised the demographic census (Fretwell & Trathan 2009, 2012, 2023):

Staining detection: emperor guano creates a dark brown stain on fast ice, visible at 30-m Landsat resolution and 10-m Sentinel-2 resolution.
Area-abundance regression: stain area correlates with colony adult count through a Fretwell calibration coefficient (~1 adult / 0.35 m\(^2\)).
High-resolution imagery: WorldView-2/3 (30 cm pixels) enables direct counting of individual birds; Maxar commercial satellites now routinely tasked during chick-rearing season.
Time-series monitoring: Sentinel-2 revisit every 5 days, allowing intra-season tracking of colony location and size.

Fretwell (2023) identified 54 colonies globally by synthesising 40 years of satellite imagery. Estimates of global emperor population total ~270 000 breeding pairs (540 000 adults), substantially higher than previous aerial-survey estimates but confirming the IUCN Vulnerable listing (2022).

\[ N_{\text{adults}} \;=\; \alpha \cdot A_{\text{stain}} + \epsilon \]

Fretwell regression: \(A_{\text{stain}}\) = guano-stain area in m\(^2\); \(\alpha \sim 2.9\) birds/m\(^2\)calibrated against validation colonies.

10. Conservation Status and Governance

Year	Designation	Authority	Scope
1959	Protected species	Antarctic Treaty	< 60°S
1964	Specially Protected Species	CCAS Protocol	< 60°S
2012	Near Threatened	IUCN Red List	Global
2022	Vulnerable	IUCN Red List	Global (uplisted)
2022	Barrows petition	US FWS	ESA listing request
2024	Threatened (ESA)	US FWS	United States jurisdiction

Ross Sea Marine Protected Area (2016)

In 2016, CCAMLR designated the Ross Sea Region Marine Protected Area — at 1.55 million km\(^2\), the world’s largest MPA. Fishery is restricted in key emperor foraging areas. A proposed MPA for the Antarctic Peninsula has been repeatedly blocked by Russia and China at CCAMLR annual meetings.

US ESA uplisting 2024

In October 2022 the US Fish and Wildlife Service proposed to list the emperor penguin as Threatened under the Endangered Species Act — the first listing of a species based primarily on projected climate impacts. The rule was finalised in 2024. Under the ESA, US agencies must consider emperor-penguin impacts in all permitting actions — including emissions-related permits, thereby establishing a novel regulatory linkage between greenhouse-gas emissions and species conservation.

11. Conservation Actions: What Can Be Done?

Unlike most conservation crises, the emperor penguin’s primary threat (climate change) cannot be addressed through local action. Feasible measures fall into three categories (Trathan et al. 2020):

1. Local stressor mitigation

Expansion of MPAs to protect foraging zones (Ross Sea MPA, proposed East Antarctica and Peninsula MPAs).
Strict biosecurity protocols on research and tourism vessels to prevent H5N1 introduction.
Limits on approach distance and visitation duration for ecotourism operators.
Precautionary catch limits on the krill fishery, with spatial closures during chick provisioning.

2. Climate-change mitigation

Global emissions pathway compatible with SSP 1-2.6 (1.5°C warming limit) is the only scenario under which emperor populations are projected to persist.
The US ESA uplisting now requires federal agencies to consider emperor impacts in greenhouse-gas permitting decisions, potentially influencing domestic emissions policy.

3. Monitoring and research

Sentinel-2 / WorldView-3 satellite census every 5 years to track colony-level demographics.
Long-term mark-recapture programmes at sentinel colonies (Pointe Géologie, Dawson-Lambton, Cape Washington).
Genomics-based sampling to quantify genetic diversity and potential climate refugia.
H5N1 surveillance and rapid-response protocols.

Crucially, translocation — a standard conservation tool for many threatened species — is not feasible for emperors. Their breeding requirements (sea ice at specific latitudes) cannot be replicated elsewhere. Captive breeding exists at a handful of aquaria but captive populations are too small and genetically narrow to function as an insurance colony.

12. Outlook: The 22nd Century Emperor

The emperor penguin of 2100 will inhabit a smaller, more fragmented range, with populations concentrated in high-latitude refugia where fast-ice persists. If current emissions trajectories continue (closer to SSP 4.5 than to SSP 2.6), the species will experience demographic collapse but likely persist in a reduced range through 2150–2200.

Critical uncertainties remain:

Colony relocation: can populations shift to suitable new sites as the climate envelope moves? Halley Bay survivors relocated ~50 km to Dawson-Lambton, but such shifts are limited by bathymetry, foraging access, and ice dynamics.
H5N1 impact: whether and when the strain reaches emperor colonies at scale is a major unknown; model mortality of 30–60% in a single outbreak is plausible.
Krill trajectory: prey-base decline interacts with direct thermal stress; compound effects could be nonlinear.
Human action: emissions scenarios remain fundamentally a policy choice. The species fate is on a timescale (decades) that is unusually amenable to adaptive management.

Every module of this course has described an adaptation tuned to a specific climatic regime: deep-diving physiology (O2 stores calibrated to a specific fast-ice-edge foraging distance), huddle thermodynamics (optimised for a specific cold-stress range), fasting biochemistry (calibrated to a specific ice-formation timetable). The emperor penguin is not merely vulnerable to climate change; the species is an indicator species for the climatic regime under which it evolved. Its persistence through 2100 will be a quantitative record of what humanity did (or did not) do about emissions.

The Groscolas 22% threshold described in Module 7 has a climatic analog: at some cumulative sea-ice-loss threshold, the species will cross a behavioural and demographic point of no return. Where that threshold sits — and whether we are on one side or the other — remains, for now, an open question.

Discussion & Graduate Exercises

Compute the population growth rate \(\lambda\) as the dominant eigenvalue of the Leslie matrix \(\mathbf{L}\) for each SSP scenario at 2050 and 2100. At what SIE does \(\lambda = 1\)?
Modify Simulation 1 to include a correlated interannual noise term (AR(1) with lag-1 coefficient 0.5 in SIE) representing ENSO/SAM teleconnection. How does autocorrelated noise affect extinction probability relative to white noise?
Using the Fretwell stain-area regression, estimate the uncertainty in global emperor population from Sentinel-2 imagery. What sample size of validation colonies is required for \(\pm 10\%\) precision?
Extend Simulation 2 to include age-structured mortality (juveniles more susceptible). Compare outbreak severity to the age-agnostic baseline.
Compute the long-term climate envelope for emperor colonies: for each 1x1° latitude-longitude grid cell around Antarctica, compute the probability that Mar–Dec fast-ice is stable for > 270 days under SSP 5-8.5 by 2100. Where are the refugia?
Estimate the carbon-emissions cost (tCO2) that would offset, through 2100, a 1-percentage-point reduction in emperor extinction probability under the Jenouvrier (2021) model.
Design a colony-level H5N1 biosecurity plan that integrates remote sensing, carcass monitoring, and human-activity restrictions. What is the cost-effectiveness (extinctions averted per dollar spent)?

Key References

• Jenouvrier, S., Holland, M., Strøeve, J., Serreze, M., Barbraud, C., Weimerskirch, H., Caswell, H. (2014). “Projected continent-wide declines of the emperor penguin under climate change.” Nature Climate Change 4, 715–718.

• Jenouvrier, S., Che-Castaño, A., Wolf, S., Holland, M., Labrousse, S., Lafond, A., Serreze, M., Sauser, C., Weimerskirch, H. (2021). “The call of the emperor penguin: legal responses to species threatened by climate change.” Global Change Biology 27, 5008–5029.

• Jenouvrier, S., et al. (2020). “The Paris Agreement objectives will likely halt future declines of emperor penguins.” Global Change Biology 26, 1170–1184.

• Fretwell, P.T., Trathan, P.N. (2019). “Emperors on thin ice: three years of breeding failure at Halley Bay.” Antarctic Science 31, 133–138.

• Fretwell, P.T., Trathan, P.N. (2009). “Penguins from space: faecal stains reveal the location of emperor penguin colonies.” Global Ecol. Biogeogr. 18, 543–552.

• Fretwell, P.T. et al. (2012). “An emperor penguin population estimate: the first global synoptic survey of a species from space.” PLOS ONE 7, e33751.

• Fretwell, P.T., Boutet, A., Ratcliffe, N. (2023). “Record low 2022 Antarctic sea ice led to catastrophic breeding failure of emperor penguins.” Communications Earth & Environment 4, 273.

• Trathan, P.N., Fretwell, P.T., Stonehouse, B. (2011). “First recorded loss of an emperor penguin colony in the recent period of Antarctic regional warming.” PLOS ONE 6, e14738.

• Trathan, P.N. et al. (2020). “The emperor penguin — vulnerable to projected rates of warming and sea ice loss.” Biol. Conserv. 241, 108216.

• Barbraud, C., Weimerskirch, H. (2001). “Emperor penguins and climate change.” Nature 411, 183–186.

• Atkinson, A., Hill, S.L., Pakhomov, E.A., et al. (2019). “Krill (Euphausia superba) distribution contracts southward during rapid regional warming.” Nat. Clim. Change 9, 142–147.

• Banyard, A.C., Lean, F.Z.X., Robinson, C., et al. (2022). “Detection of highly pathogenic avian influenza virus H5N1 clade 2.3.4.4b in great skuas.” Viruses 14, 212.

• Ashcroft, R., et al. (2024). “Spread of H5N1 highly pathogenic avian influenza in Antarctic seabirds.” Nat. Commun. 15, 7042.

• Bustamante, P., et al. (2021). “Mercury contamination in Antarctic emperor penguins over a century.” Environ. Pollut. 281, 117026.

• Le Guen, C., et al. (2018). “Foraging behaviour and annual oceanic provisioning strategy of emperor penguins.” Mar. Ecol. Prog. Ser. 597, 191–206.

• Parkinson, C.L. (2019). “A 40-year record reveals gradual Antarctic sea ice increases followed by decreases.” PNAS 116, 14414–14423.

• US Fish and Wildlife Service (2024). “Endangered and Threatened Wildlife and Plants; Threatened Species Status for the Emperor Penguin.” Federal Register 89 FR 83834.

• Caswell, H. (2001). Matrix Population Models: Construction, Analysis, and Interpretation, 2nd ed. Sinauer.

• IUCN (2022). “Aptenodytes forsteri: emperor penguin uplisted to Vulnerable.” IUCN Red List assessment.

Course Synthesis

Across eight modules we have traced the emperor penguin from evolutionary origin (Module 0) through thermoregulation (Module 1), huddle dynamics (Module 2), breeding physiology (Module 3), deep-diving oxygen management (Module 4), underwater hydrodynamics (Module 5), two-voice vocal recognition (Module 6), and the 120-day fasting biochemistry (Module 7). In every system, the bird operates at the edge of what a 30-kg vertebrate can plausibly do. Each adaptation is tuned to an environment — the seasonal fast-ice cycle, the chorus of a 10 000-bird colony, the dark pressure of 500-m dives — that is changing.

The question this final module leaves open is whether the emperor penguin, a species defined by its relationship with Antarctic sea ice, will survive the rapid decoupling of that relationship in the 21st century. The demographic models say: only under SSP 1-2.6. Every other trajectory implies steep decline. The tool to deliver SSP 1-2.6 lies, uniquely among conservation problems, outside the Antarctic: in the global decisions about emissions.

We close where Le Maho began in 1977: the emperor penguin is a strategy to live and breed in the cold. Whether the cold remains is now up to us.

← Module 7: Fasting Biochemistry Back to Course Overview →