Setting up the environment

The first step is to set up the conda environment. You should use the profsea-env.yml environment file, building it with

conda env create -f profsea-env.yml

and then activating with conda activate profsea.

---

One final step here is to install ProFSea as an editable package, using

pip install -e .

while inside the highest level directory.

Importing components

The next step is to import all the model components you want to run simulations with - a full list is available in the documention `here <>`__. This can be any combination of available global and spatial components, how you build your model is up to you!

To generate a comprehensive set of sea-level projections, we’ll run with all the components that contribute to sea-level change:

• thermal expansion

• glaciers

• greenland ice-sheet

• antarctic ice-sheet

• landwater

In the example below, we’ll run with the configuration being used for SLEIP and Munday et al. (2026), in prep.

from profsea.components.core.global_model import Global
from profsea.components.global_ import (
    LandwaterAR6,
    GreenlandAR6,
    ThermalExpansion,
    AntarcticaISMIP6,
    Glacier,
)

ProFSea can run with a single or an ensemble of temperature and ocean heat content forcing trajectories. The trajectories need to be anomalies, relative to some baseline period - ProFSea assumes by default this is 1996-2014

Here we’ll just run with an idealised temperature ramp up and stabilisation, and corresponding ocean heat content change, running from 2006 to 2300:

import matplotlib.pyplot as plt
import numpy as np

t_change = np.concatenate(
    [np.linspace(0.01, 8, 150), np.linspace(8, 8, 145)]
)
t_change = t_change[None, :] * np.random.normal(loc=1.0, scale=0.1, size=(100, 1)) # add some variability across members
ohc_change = t_change * 1e24 # in Joules
years = np.arange(2006, 2301)
years = np.broadcast_to(years, t_change.shape)

plt.figure(figsize=(4, 2))
plt.plot(years, t_change, color="royalblue", alpha=0.1)
plt.plot(years.mean(axis=0), t_change.mean(axis=0), color="royalblue", label="GMST change")
plt.xlabel("Year")
plt.ylabel("GMST change (°C)")
plt.show()
plt.close()

Building the model

To build your model, simply add your components to a dictionary, with whatever names you like, and pass it to the Global() class

global_components = {
    "landwater": LandwaterAR6(),
    "greenland": GreenlandAR6(),
    "expansion": ThermalExpansion(ohc_change), # only the thermal expansion needs OHC change, so we'll pass it in here
    "wais": AntarcticaISMIP6(ais_region="wais"),
    "eais": AntarcticaISMIP6(ais_region="eais"),
    "peninsula": AntarcticaISMIP6(ais_region="peninsula"),
    "glacier": Glacier(),
}
global_model = Global(components=global_components, end_yr=2301, num_members=1000)

now we’re ready to run our simulation…

projections = global_model.run(
    T_change=t_change,
    scenario="stabilisation",
    member_seed=42
)
global_model.sum_components(projections)
gmslr = global_model.results["gmslr"]

visualise the global simulations!

plt.figure(figsize=(4, 3))
lower = np.percentile(gmslr, 5, axis=0)
upper = np.percentile(gmslr, 95, axis=0)
plt.fill_between(years[0], lower, upper, color="lightcoral", alpha=0.5, label="90% range")
plt.plot(years[0], gmslr.mean(axis=0), color="darkred")
plt.xlabel("Year")
plt.ylabel("GMSLR (m)")
plt.show()
plt.close()

Building the Spatial model

Start with the imports:

from profsea.components.core import Spatial
from profsea.components.spatial import (
    SterodynamicCMIP6,
    Fingerprint,
    GIA,
)

Now we’ll build the model in the same way as the global model, but we’ll pass the individual components’ global projections to each spatial module.

spatial_components = {
    "sterodynamic": SterodynamicCMIP6(
        projections["expansion"], # passing in the global thermal expansion projection
    ),
    "greenland": Fingerprint(
        projections["greenland"], fingerprint_component="greenland"
    ),
    "landwater": Fingerprint(
        projections["landwater"],
        fingerprint_component="landwater",
    ),
    "wais": Fingerprint(
        projections["wais"],
        fingerprint_component="wais",
    ),
    "eais": Fingerprint(
        projections["eais"],
        fingerprint_component="eais",
    ),
    "glacier": Fingerprint(
        projections["glacier"],
        fingerprint_component="glacier",
    ),
    "gia": GIA(
        sample_spatial=False,
    ),
}

# Pass to the spatial model
spatial_model = Spatial(components=spatial_components)
spatial_model.run(member_seed=42)

total_rsl = spatial_model.sum_components(spatial_model.results)

[13:46:57] ✓ ProFSea assets found locally!                                                     spatial_model.py:359

           Baseline period = 1995 to 2014                                                      spatial_model.py:107

           Simulating 7 sea-level components...: sterodynamic, greenland, landwater, wais,     spatial_model.py:188
           eais, glacier, gia                                                                                      

Saving the spatial projections is easy + quick using Zarr format:

spatial_model.save_components(
    spatial_model.results, scenario_name="stabilisation", output_format="zarr"
)

[13:47:42] ✓ Successfully saved Zarr: ./stabilisation_spatial_projection.zarr                  spatial_model.py:331

           Output shape was (5, 295, 180, 360) (members, time, lat, lon)                       spatial_model.py:335

Plot the output! This takes about 90 seconds, since our data goes from lazy → in memory

import cartopy.crs as ccrs

fig = plt.figure(figsize=(6, 4), layout="constrained")

total_rsl_med = total_rsl[2, -1, :, :]
vmax = np.nanmax(np.abs(total_rsl))
vmin = -vmax

ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
ax.set_title("50th Percentile, Year 2300")
ax.pcolormesh(
    spatial_model.grid_lons,
    spatial_model.grid_lats,
    total_rsl_med,
    transform=ccrs.PlateCarree(),
    cmap="PuOr_r",
    vmin=vmin,
    vmax=vmax,
)
ax.coastlines()
fig.colorbar(
    mappable=ax.collections[0],
    label="Relative Sea Level (m)",
    orientation="horizontal",
    pad=0.02,
)
plt.show()