Sea level rise data explained

On this page

Problem statement

"Sea level rise" is at least four different quantities, and a map that mixes them is wrong in ways that no validation catches. Global mean sea level rise, regional sea level change, relative sea level change at a gauge and the local flood level in a given year are related and not interchangeable, and the differences are large enough to change which streets appear on a map.

The commonest error is simpler than any of that: adding a projected rise to a terrain model that is in a different vertical datum. A national orthometric datum, a tidal datum and the ellipsoid are tens of centimetres to tens of metres apart, and a rise of half a metre added to the wrong zero produces a map of something else entirely.

Quick answer

Get every component into one vertical datum before adding anything:

FLOOD_LEVEL_M = (
    mhhw_in_navd88          # the tidal surface, converted through the gauge datums
    + projected_rise_m      # the scenario, relative to a stated baseline period
    + storm_allowance_m     # the return-period surge, if you are mapping a flood
    - vertical_land_motion  # subsidence adds, uplift subtracts — unless the
)                           # projection is already relative, as NOAA's station ones are

Every term needs a datum, a baseline period and a source. At Boston, MHHW is 1.453 m above NAVD88, and NOAA's 2022 Intermediate scenario for the gauge is +1.13 m by 2100, so the still-water level is about 2.58 m NAVD88 before any storm allowance (Example 1).

The terms that make up a mapped coastal flood level. assembling a flood leveltidal datumMHHW 1.453 m NAVD88local, 19-year epoch+ projected risescenario and yearrelative to a baseline+ storm allowancereturn periodwhat actually floods− vertical land motionmm/yr × yearssubsidence adds= level in the DEM datumone numbernow it can be mappedFour terms, four datums, four sources spatialworkflow.io
Four terms, four datums and four sources; the sum is only as good as the worst of them.

Four quantities, land motion, baselines and scenarios

1. Distinguish the four quantities

  • Global mean sea level — the ocean-volume average from satellite altimetry since 1993 and from tide gauges before that. A single global number.
  • Regional sea level change — the same quantity with ocean dynamics, gravity and ice-mass fingerprints included. Varies by a factor of two or more between regions.
  • Relative sea level — what a tide gauge measures: the sea relative to the land it is bolted to, so it includes vertical land motion.
  • Extreme water level — the still water plus the surge and wave setup for a given return period. This is what floods anything.

A projection published as "global mean" is not the number to add at a specific place.

2. Get the vertical land motion right

Land subsidence adds directly to relative sea level rise and can exceed the climatic signal. Subsiding deltas and formerly glaciated coasts move in opposite directions, and GNSS or InSAR rates are the source.

3. Read the baseline period of the projection

A projection is relative to a baseline year, and a tidal datum to its epoch. NOAA's 2022 scenarios for Boston start from 2005; Boston's tidal datums are for the 1983–2001 epoch, centred on 1992. Adding one to the other leaves out the rise between the two — thirteen years at the gauge's 2.97 mm a year, about 4 cm. The error runs the other way when a projection's baseline is earlier than the datum's epoch: then the overlap is counted twice. Either way, line the two up before adding.

4. Choose a scenario and name it

Projections come as scenarios with likely ranges and low-confidence high-end branches. A map showing one number without naming the scenario, the percentile and the year is not interpretable.

5. Convert the tidal surface into the terrain model's datum

This is the step that gets skipped. Terrain models are in a national orthometric datum; tidal datums are local and published relative to station datum. Convert through the gauge's published offsets.

6. Add a storm allowance if you are mapping flooding

Still-water sea level rise alone floods almost nothing new. The damaging events are surges on top of a high tide, and the return-period level is what coastal flood maps are actually drawn at.

7. Publish the uncertainty

Every term has one: the projection range, the land motion rate, the DTM's vertical error, and the datum conversion. A single contour implies a precision none of them supports.

Global mean, regional, relative and extreme water level compared. four different quantitiesglobal meanrelativeextreme levelocean volume✓yes✓yes✓yesregional pattern✕no✓yes✓yesland motion✕no✓yes✓yessurge and waves✕no✕no✓yesOnly the last one floods anything spatialworkflow.io
Four quantities, four uses; only the last one floods anything.

A Boston flood level with every term named

The data used on this page

All for Boston: the harbour DTM and tidal datums from extracting a coastline from a raster, and three of NOAA's derived products for the Boston gauge — its long-term relative sea level trend, the 2022 interagency sea level rise scenarios for the station, and its extreme water level statistics.

import numpy as np
import rasterio
import requests
from rasterio.enums import Resampling
from rasterio.warp import calculate_default_transform, reproject, transform_bounds

BBOX = (-71.08, 42.30, -70.93, 42.40)                     # Boston Harbor, lon/lat
SRC = "https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/current/n43w072/USGS_13_n43w072.tif"

with rasterio.open(SRC) as src:
    window = rasterio.windows.from_bounds(*transform_bounds("EPSG:4326", src.crs, *BBOX),
                                          transform=src.transform).round_offsets().round_lengths()
    data = src.read(1, window=window)
    src_transform, src_crs = src.window_transform(window), src.crs
dst_transform, width, height = calculate_default_transform(
    src_crs, "EPSG:26919", data.shape[1], data.shape[0],
    *rasterio.transform.array_bounds(*data.shape, src_transform), resolution=10)
dtm = np.full((height, width), np.nan, dtype="float32")
reproject(data, dtm, src_transform=src_transform, src_crs=src_crs, dst_transform=dst_transform,
          dst_crs="EPSG:26919", resampling=Resampling.bilinear, src_nodata=-999999, dst_nodata=np.nan)
with rasterio.open("dtm_navd88.tif", "w", driver="GTiff", width=width, height=height, count=1,
                   dtype="float32", crs="EPSG:26919", transform=dst_transform, nodata=np.nan) as out:
    out.write(dtm, 1)

datums = requests.get("https://api.tidesandcurrents.noaa.gov/mdapi/prod/webapi/stations/8443970/"
                      "datums.json", params={"units": "metric"}, timeout=60).json()
station = {d["name"]: d["value"] for d in datums["datums"]}
NAVD88 = {name: round(station[name] - station["NAVD88"], 3) for name in ("MHHW", "MSL", "MLLW")}
print(f"DTM {width} x {height} cells of 10 m, {np.nanmin(dtm):.1f} to {np.nanmax(dtm):.1f} m NAVD88")
print(f"Boston tidal datums in NAVD88 (epoch {datums['epoch']}): {NAVD88}")

# NOAA's derived products for the same gauge: trend, 2022 scenarios, extreme water levels
DP = "https://api.tidesandcurrents.noaa.gov/dpapi/prod/webapi"
trend = requests.get(f"{DP}/product/sealvltrends.json", params={"station": "8443970"}, timeout=60).json()["SeaLvlTrends"][0]
TREND_MM_YR = trend["trend"] * 25.4 / 10                      # published in inches a decade
projections = requests.get(f"{DP}/product/slr_projections.json", params={"station": "8443970"},
                           timeout=60).json()["SlrProjections"]
RISE_2100_M = {p["scenario"]: p["projectionRsl"] / 100 for p in projections if p["projectionYear"] == 2100}
BASELINE_YEAR = min(p["projectionYear"] for p in projections)
exceedance = requests.get(f"{DP}/product.json", timeout=60, params={
    "name": "extremewaterlevels", "station": "8443970", "details": "exceedancelevels",
    "levelType": "high", "units": "metric"}).json()
ONE_PERCENT_ABOVE_MHHW = next(e["level"] for e in exceedance["GEV_HI"]["exceedanceLevels"]
                              if e["probability"] == 0.01)
print(f"relative trend {TREND_MM_YR:.2f} mm/yr ({trend['startDate']}–{trend['endDate']})")
print(f"2022 scenarios for 2100, relative to {BASELINE_YEAR}: {RISE_2100_M}")
print(f"1% annual-chance level: {ONE_PERCENT_ABOVE_MHHW} m above MHHW")
DTM 691 x 1127 cells of 10 m, -9.8 to 49.3 m NAVD88
Boston tidal datums in NAVD88 (epoch 1983-2001): {'MHHW': 1.453, 'MSL': -0.092, 'MLLW': -1.678}
relative trend 2.97 mm/yr (01/15/1921–12/15/2025)
2022 scenarios for 2100, relative to 2005: {'Low': 0.48, 'Intermediate-Low': 0.67, 'Intermediate': 1.13, 'Intermediate-High': 1.48, 'High': 1.9}
1% annual-chance level: 1.41 m above MHHW

Example 1 — assemble the flood level with every term named

import dataclasses, json

@dataclasses.dataclass
class FloodLevel:
    tidal_datum: str              # e.g. "MHHW"
    tidal_datum_m_navd88: float   # from the gauge's published offsets
    scenario: str                 # e.g. "intermediate, 50th percentile"
    baseline_period: str          # e.g. "2000"
    projected_rise_m: float
    year: int
    storm_return_period_years: int | None
    storm_allowance_m: float
    vertical_land_motion_mm_yr: float
    reference_year: int = 2026

    def level_m_navd88(self):
        vlm = (self.vertical_land_motion_mm_yr / 1000) * (self.year - self.reference_year)
        return (self.tidal_datum_m_navd88 + self.projected_rise_m
                + self.storm_allowance_m - vlm)

f = FloodLevel(tidal_datum="MHHW", tidal_datum_m_navd88=NAVD88["MHHW"],
               scenario="NOAA 2022 Intermediate, relative sea level", baseline_period=str(BASELINE_YEAR),
               projected_rise_m=RISE_2100_M["Intermediate"], year=2100,
               storm_return_period_years=100, storm_allowance_m=ONE_PERCENT_ABOVE_MHHW,
               vertical_land_motion_mm_yr=0.0)    # already inside a relative projection
print(json.dumps(dataclasses.asdict(f) | {"level_m_navd88": round(f.level_m_navd88(), 3)},
                 indent=2))
{
  "tidal_datum": "MHHW",
  "tidal_datum_m_navd88": 1.453,
  "scenario": "NOAA 2022 Intermediate, relative sea level",
  "baseline_period": "2005",
  "projected_rise_m": 1.13,
  "year": 2100,
  "storm_return_period_years": 100,
  "storm_allowance_m": 1.41,
  "vertical_land_motion_mm_yr": 0.0,
  "reference_year": 2026,
  "level_m_navd88": 3.993
}

A dataclass rather than a float, so the number cannot travel without its assumptions. The land-motion term is zero here on purpose: NOAA's station scenarios are relative sea level, with the gauge's background rate already inside them, and subtracting land motion again would count it twice. Use the term only with a global or absolute projection — and then a subsiding coast has a negative rate, because land going down raises relative sea level.

Example 2 — the trend from a real gauge record

import urllib.request, json, numpy as np, pandas as pd

def monthly_means(station, begin_year, end_year):
    frames = []
    for y in range(begin_year, end_year + 1):
        url = (f"https://api.tidesandcurrents.noaa.gov/api/prod/datagetter"
               f"?product=monthly_mean&application=spatialworkflow-docs"
               f"&begin_date={y}0101&end_date={y}1231&datum=MSL&station={station}"
               f"&time_zone=gmt&units=metric&format=json")
        with urllib.request.urlopen(url, timeout=120) as r:
            raw = json.loads(r.read())
        if "data" in raw:
            frames.append(pd.DataFrame(raw["data"]))
    df = pd.concat(frames, ignore_index=True)
    df["t"] = pd.to_datetime(df["year"] + "-" + df["month"].str.zfill(2) + "-01")
    df["msl"] = pd.to_numeric(df["MSL"], errors="coerce")
    return df.dropna(subset=["msl"])

df = monthly_means("8443970", 1990, 2025)
x = (df["t"] - df["t"].min()).dt.days / 365.25
slope, _ = np.polyfit(x, df["msl"], 1)
print(f"{len(df)} monthly means, trend {slope * 1000:.2f} mm/yr")
print(f"NOAA's published trend since {trend['startDate'][-4:]}: {TREND_MM_YR:.2f} mm/yr")
429 monthly means, trend 4.98 mm/yr
NOAA's published trend since 1921: 2.97 mm/yr

The two numbers differ because the rate is not constant: since 1990 the gauge has risen at 4.98 mm a year against 2.97 over the full record from 1921. A trend is a statement about a period, and has to be quoted with it.

A gauge trend is relative sea level: it includes whatever the land is doing. Comparing it with a global mean number without separating vertical land motion is the classic mistake.

Example 3 — a scenario sweep rather than a single line

import numpy as np, rasterio
from rasterio.features import shapes
from shapely.geometry import shape

SCENARIOS = RISE_2100_M                  # NOAA 2022, Boston, 2100

with rasterio.open("dtm_navd88.tif") as src:
    dtm = src.read(1, masked=True)
    transform, crs = src.transform, src.crs

base = NAVD88["MHHW"]
cell_km2 = abs(transform.a * transform.e) / 1e6
today = (dtm <= base).filled(False)
for name, rise in SCENARIOS.items():
    level = base + rise
    below = (dtm <= level).filled(False)
    print(f"{name:18} +{rise:.2f} m → level {level:.3f} m NAVD88, "
          f"{(below & ~today).sum() * cell_km2:5.1f} km² of land newly below MHHW")
Low                +0.48 m → level 1.933 m NAVD88,   1.1 km² of land newly below MHHW
Intermediate-Low   +0.67 m → level 2.123 m NAVD88,   1.7 km² of land newly below MHHW
Intermediate       +1.13 m → level 2.583 m NAVD88,   4.7 km² of land newly below MHHW
Intermediate-High  +1.48 m → level 2.933 m NAVD88,   9.6 km² of land newly below MHHW
High               +1.90 m → level 3.353 m NAVD88,  17.0 km² of land newly below MHHW

This is a bathtub count: every 10 m cell below the level, whether or not water could reach it. It is an upper bound for the still-water case and an underestimate of flooding, which comes with the 1% storm on top.

Publish the sweep, not one scenario. The difference between the low and high lines is the honest message, and a single contour hides it.

Land in the Boston Harbor DTM newly below mean higher high water in 2100 under the five NOAA 2022 scenarios for the Boston gauge. Boston, 2100: five scenarios, one DTMLow, +0.48 m1.1 km²Int.-Low, +0.67 m1.7 km²Intermediate, +1.13 m4.7 km²Int.-High, +1.48 m9.6 km²High, +1.90 m17.0 km²09 km²18 km²A fifteen-fold range from one gauge, one year and one terrain model spatialworkflow.io
Measured in Example 3, a bathtub count on the 3DEP DTM. The spread between scenarios is the message.

Why global and local numbers differ

Why global and local numbers differ so much

Sea level is not flat. Ocean currents, thermal expansion patterns, and the gravitational fingerprint of melting ice sheets — meltwater from Greenland raises sea level less nearby and more in the far field — all mean regional change departs substantially from the global mean. Add vertical land motion and the relative change at a gauge can be several times the global figure or, on rebounding coasts, negative.

Why still-water rise alone maps almost nothing

Half a metre of still-water rise inundates only land below half a metre above the present high-water line, which is a narrow band. What changes is the frequency of extremes: an event that used to be a 1-in-100-year level becomes a 1-in-10 or annual one. That is the honest way to present the result, and it is why a coastal flood map is drawn at a return-period level rather than at the rise itself.

Why the baseline period is a real source of double counting

Both the tidal datum and the projection are anchored to a period. A datum from the 1983–2001 epoch already reflects the sea level of the 1990s; a projection relative to 2000 counts from roughly the same point, so adding them is approximately right. A projection relative to 1900 added to a modern datum counts the twentieth century twice.

Why the terrain model's error belongs in the total

A LiDAR DTM has a vertical RMSE of 0.1–0.3 m in open terrain and worse in vegetation and on soft ground. Against the Low scenario's +0.48 m for Boston that is a large fraction of the signal, which is why inundation extents should be presented as a band rather than a line.

What still-water rise does against what it does to the frequency of extremes. the honest framingstill-water rise alone✓inundates a narrow band✓below the new high-water line✓looks unimpressive on a map✓is not what floods anythingchange in frequency✕a 1-in-100 level becomes 1-in-10✕then annual✕same level, different return period✕what flood policy is written inMap return-period levels, not the rise itself spatialworkflow.io
The damaging change is in the frequency of extremes, not in the mean.

Relative versus absolute, and altimetry's start

  • Relative and absolute differ by the land motion. Say which you are using.
  • Tide gauge records need datum continuity. Gauges are moved and re-levelled.
  • Altimetry starts in 1993. Earlier global series are gauge reconstructions.
  • Scenario percentiles matter. A 50th and a 95th differ by a factor of two or more.
  • Wave setup and runup are on top of surge, and are not in a still-water level.
  • Groundwater and drainage flood inland before the sea reaches there.
  • Defences are not in a DTM. A bathtub map floods behind seawalls.
  • State every datum. The commonest error is arithmetic between different zeros.

FAQ

What is the difference between global and relative sea level rise?

Global mean is the ocean-volume average. Relative is what a gauge measures at a place, including regional ocean dynamics and vertical land motion, and it can be several times the global figure.

Can I add a projected rise to my DEM?

Only after converting both to the same vertical datum, and only after adding the tidal surface. A rise added to a bare DEM zero is a map of nothing in particular.

Why does still-water rise flood so little?

Because it only inundates the narrow band below the new high-water line. The real effect is that extreme levels become far more frequent, which is why flood maps use return-period levels.

Which scenario should I map?

Several, and name each one with its percentile and year. A single contour implies a certainty that no projection supports.

Does land subsidence matter?

Often more than the climatic signal. In subsiding deltas it dominates relative sea level rise; on rebounding coasts it partly cancels it.

How accurate is the resulting extent?

No better than the DTM, which is typically 0.1–0.3 m RMSE — a large fraction of the Low scenario's half metre. Present a band, not a line.