Oil slicks near the Santa Ynez Unit platforms before and after restart

A before/after/control/impact (BACI) look at SkyTruth Cerulean detections in the Santa Barbara Channel

NoteBottom line

After Harmony’s test restart (May 1, 2025), the share of sea surface with a detected slick within 5 km of the three platforms rose from 0.56% to 1.53% of imaged area per satellite pass. Over the rest of the Channel it rose from 0.72% to 0.89%. The increase is concentrated at Harmony and Hondo; Heritage shows no change.

The signal is suggestive but not conclusive. After adjusting for season, the effect is +0.66 percentage points (95% CI -0.07 to 1.43). About 8% of randomly placed “placebo” platform zones in the Channel show a change at least as large. The result also depends on the comparison area: against the nearby Coal Oil Point seep band, which also got oilier, the effect is close to zero.

1 Question and background

Did the restart of oil production at platforms Harmony, Heritage and Hondo (Santa Ynez Unit, operated by Sable Offshore) change the number or area of oil slicks detected in the Santa Barbara Channel?

The platforms were shut in after the 2015 Refugio pipeline spill. The restart was staggered:

Restart timeline (see Sources)

Date	Event
May 2025	Sable restarts six wells on Harmony as a test. This analysis uses 2025-05-01 as the start of the “After” period.
2026-03-14	Shipping from Las Flores Canyon resumes under a Department of Energy order. Harmony produces about 22,000 b/d. This starts phase 2 here.
~April 2026	Heritage cleared by BSEE to resume operations.
Q3 2026 (expected)	Hondo expected to come online.

Because of the staggered timeline, results are shown for the whole After period and for two phases: phase 1 (Harmony test only) and phase 2 (full restart).

2 Data

Slicks. All slick polygons from SkyTruth Cerulean’s public API (public.slick_plus) for the study area and dates in the settings:

https://api.cerulean.skytruth.org/collections/public.slick_plus/items?f=csv&limit=10000&offset=0&datetime=2023-01-01T00:00:00Z/2026-12-31T23:59:59Z&bbox=-120.65,33.85,-119,34.5

That returned 1,235 slicks from Jan 2023 to Sep 2026.

Imaging effort. Cerulean only reports scenes where it found a slick, but a slick can only be seen when a Sentinel-1 radar satellite passes overhead. I therefore also downloaded every Sentinel-1 scene Cerulean ingested over the same box (public.sentinel1_grd), including scenes with no detections. Adjacent frames from one orbit were merged into a single pass. That gives 791 passes, 312 of which image at least 90% of the platform zone and are used in the BACI comparisons.

# Download slicks and satellite scenes from the Cerulean API, unless we already
# have a saved copy in data/ and `refresh_data` is FALSE.
bbox_str  <- paste(bbox[c("xmin", "ymin", "xmax", "ymax")], collapse = ",")
dates_str <- sprintf("%sT00:00:00Z/%sT23:59:59Z", date_start, date_end)
api_url <- function(collection, format, limit, offset)
  sprintf("%s/%s/items?f=%s&limit=%d&offset=%d&datetime=%s&bbox=%s",
          cerulean_api, collection, format, limit, offset, dates_str, bbox_str)

if (refresh_data || !file.exists("data/slicks.csv") || !file.exists("data/s1_scenes.json")) {
  dir.create("data", showWarnings = FALSE)

  # Slicks (public.slick_plus): CSV, fetched 10,000 rows at a time
  pages <- list(); offset <- 0
  repeat {
    tmp <- tempfile(fileext = ".csv")
    download.file(api_url("public.slick_plus", "csv", 10000, offset), tmp, quiet = TRUE)
    page <- read_csv(tmp, col_types = cols(.default = "c"))
    pages[[length(pages) + 1]] <- page
    if (nrow(page) < 10000) break
    offset <- offset + 10000
  }
  bind_rows(pages) %>% distinct(id, .keep_all = TRUE) %>% write_csv("data/slicks.csv", na = "")

  # Every Sentinel-1 scene Cerulean processed over the box (public.sentinel1_grd),
  # including scenes with no slicks: this is our measure of imaging effort.
  pages <- list(); offset <- 0
  repeat {
    page <- fromJSON(api_url("public.sentinel1_grd", "json", 500, offset))
    pages[[length(pages) + 1]] <- page
    if (NROW(page) < 500) break
    offset <- offset + 500
  }
  write_json(bind_rows(pages), "data/s1_scenes.json")
}

# Divide the ocean in the study area into square cells (default 1 km).
crs_m <- 3310  # California Albers, in metres

platforms <- platforms_ll %>% st_as_sf(coords = c("lon", "lat"), crs = 4326) %>% st_transform(crs_m)

# Land (Natural Earth 1:50m); ocean = study box minus land
land <- rnaturalearth::ne_countries(scale = 50, returnclass = "sf") %>% st_make_valid() %>%
  st_crop(st_bbox(bbox + c(-1, -1, 1, 1), crs = 4326)) %>% st_union() %>% st_transform(crs_m)
domain <- st_as_sfc(st_bbox(bbox, crs = 4326)) %>% st_transform(crs_m)
ocean  <- st_difference(domain, land)

# Keep cells whose centre is in the ocean
grid <- st_make_grid(domain, cellsize = cell_m) %>% st_sf(geometry = .)
grid <- grid[lengths(st_intersects(st_centroid(grid), ocean)) > 0, ]
grid$cell    <- seq_len(nrow(grid))
grid$cell_m2 <- as.numeric(st_area(grid))
centres <- st_centroid(grid)

# Distance to the nearest platform, and which platform that is
dist_to_platforms <- st_distance(centres, platforms)
grid$dist_km <- as.numeric(apply(dist_to_platforms, 1, min)) / 1000
grid$nearest <- platforms$name[apply(dist_to_platforms, 1, which.min)]

# Cell coordinates (metres and longitude) and distance to the mainland coast
xy <- st_coordinates(centres); grid$x <- xy[, 1]; grid$y <- xy[, 2]
grid$lon <- st_coordinates(st_transform(centres, 4326))[, 1]
land_parts <- st_cast(land, "POLYGON")
mainland   <- land_parts[which.max(st_area(land_parts))]
grid$coast_km <- as.numeric(st_distance(centres, mainland)) / 1000

# Build the core dataset: for every satellite pass, which cells it imaged and how
# much of each cell was oiled. Cells that were imaged but had no slick get oil = 0,
# which is what lets us compare rates rather than raw counts.

# ---- Satellite scenes and their footprints
scenes <- fromJSON("data/s1_scenes.json") %>% as_tibble() %>%
  distinct(scene_id, .keep_all = TRUE) %>%
  mutate(t = ymd_hms(start_time), sat = substr(scene_id, 1, 3)) %>%   # sat = S1A, S1C, S1D
  filter(t >= as.POSIXct(date_start, tz = "UTC"), t < as.POSIXct(date_end, tz = "UTC") + 86400)
scene_geom <- st_as_sfc(sub("^SRID=4326;", "", scenes$geometry), crs = 4326) %>% st_transform(crs_m)
scenes <- st_sf(scenes %>% select(-geometry, -url), geometry = scene_geom)
scenes$scene <- seq_len(nrow(scenes))

# A cell counts as imaged by a scene when the cell's centre is inside the scene footprint
covered  <- st_intersects(scenes, centres)
scenes   <- scenes[lengths(covered) > 0, ]
covered  <- covered[lengths(covered) > 0]
coverage <- data.table(scene = rep(scenes$scene, lengths(covered)), cell = unlist(covered))

# ---- Slicks, matched to the scene they were detected in
slicks <- read_csv("data/slicks.csv", show_col_types = FALSE)
n_slicks_downloaded <- nrow(slicks)
slick_geom <- st_as_sfc(sub("^SRID=4326;", "", slicks$geometry), crs = 4326) %>% st_make_valid() %>% st_transform(crs_m)
slicks <- st_sf(slicks %>% select(id, slick_timestamp, machine_confidence, cls, hitl_cls, area, length,
                                  s1_scene_id, source_type_2_ids),
                geometry = slick_geom) %>%
  mutate(t = ymd_hms(slick_timestamp))
slicks$scene <- scenes$scene[match(slicks$s1_scene_id, scenes$scene_id)]
slicks <- slicks[!is.na(slicks$scene), ]

# Oiled area of each slick inside each grid cell
pieces <- st_intersection(slicks %>% select(id, scene), grid %>% select(cell))
slick_cells <- data.table(id = pieces$id, scene = pieces$scene, cell = pieces$cell,
                          oil_m2 = as.numeric(st_area(pieces)))
oil <- slick_cells[, .(oil_m2 = sum(oil_m2)), by = .(scene, cell)]

# Scene x cell panel: every imaged cell, plus oil in cells whose centre falls just
# outside the footprint edge (the slick proves the cell was seen)
scene_cells <- merge(coverage, oil, by = c("scene", "cell"), all.x = TRUE)
scene_cells[is.na(oil_m2), oil_m2 := 0]
scene_cells <- rbind(scene_cells, oil[!coverage, on = .(scene, cell)])

# ---- Merge adjacent frames of one orbit into a single satellite "pass"
scene_dt <- as.data.table(st_drop_geometry(scenes))
scene_dt[, pass := paste(sat, absolute_orbit_number)]
scene_dt[, t_pass := min(t), by = pass]
passes <- unique(scene_dt[, .(pass, t = t_pass, sat)])
passes[, period := fcase(t < restart, "Before", t < restart2, "After: phase 1", default = "After: phase 2")]
passes[, month := floor_date(t, "month")]

# cell_obs: one row per pass x imaged cell. Frames overlap slightly, so take the max.
cell_obs <- merge(scene_cells, scene_dt[, .(scene, pass)], by = "scene")[
  , .(oil_m2 = max(oil_m2)), by = .(pass, cell)]
cell_obs[, cell_m2 := grid$cell_m2[cell]]

# Slick attributes plus the pass and period each slick belongs to
slick_dt <- as.data.table(st_drop_geometry(slicks)) %>%
  merge(scene_dt[, .(scene, pass)], by = "scene") %>%
  merge(passes[, .(pass, period)], by = "pass")

Table 1: Sentinel-1 passes over the study area by year and satellite. Sentinel-1B failed in Dec 2021; Sentinel-1C entered service in spring 2025 and Sentinel-1D in early 2026.

Year	S1A	S1C	S1D	Total
2023	147	0	0	147
2024	148	0	0	148
2025	148	114	0	262
2026	72	96	66	234

WarningThe main confounder: satellite coverage doubled at the same time as the restart

Sentinel-1C came online within weeks of the Harmony restart. Raw slick counts therefore roughly double everywhere in the Channel, whatever the platforms did (Figure 1). Every comparison below is per satellite pass, and compares the platform zone against a control area imaged on the same passes.

Figure 1: Top: satellite passes per quarter. Bottom: raw slick detections per quarter. Dashed lines mark the Harmony test restart and the full restart.

3 Methods

1. Grid. The ocean part of the study area is divided into 1 km cells (California Albers projection; land from Natural Earth 1:50m).
2. Panel. For every pass, I record which cells it imaged (cell centre inside the scene footprint) and the oiled area of every slick in each cell. A cell counts as imaged on a pass even when no slick was found, which puts the zeros back into the data.
3. Zones. Impact = cells within 5 km of any of the three platforms. Control = cells more than 10 km from every platform. Cells in between are dropped to reduce spillover.
4. Outcomes per pass and zone.
   • Occurrence: share of imaged cells touched by at least one slick.
   • Oiled fraction: oiled area ÷ imaged area.
   Both are rates per pass, so they don’t depend on how often the Channel is imaged.
5. BACI estimate. Difference-in-differences: (impact After − impact Before) − (control After − control Before). A ratio-of-ratios version is also shown.
6. Uncertainty and robustness.
   • Seasonal adjustment. The After period over-represents May–September, so pass-level (impact − control) differences are regressed on a post indicator plus month-of-year fixed effects.
   • Block bootstrap. 95% intervals come from resampling calendar months within each period, which keeps weather-correlated passes together.
   • Placebo zones. The three-platform layout is shifted across the Channel, and the effect is recomputed at every offset that stays over water and doesn’t overlap the real platforms. This shows how unusual the real change is.
   • Placebo restart date. A fake restart one year earlier, using only pre-restart data.
   • Sensitivity checks: zone radius, choice of control area, and excluding slicks Cerulean classifies as natural seeps.

# The building blocks every estimate below uses.

# Label cells as "Impact", "Control" or NA (dropped) by distance to the nearest platform
zone_of <- function(dist_km, r_imp = R_imp, r_excl = R_excl)
  fifelse(dist_km <= r_imp, "Impact", fifelse(dist_km > r_excl, "Control", NA_character_))
zone_map <- function(zone) data.table(cell = grid$cell, zone = zone)[!is.na(zone)]

# For each pass: occurrence and oiled fraction in the impact and control zones.
# Returns one row per usable pass (impact zone >= min_imp_cov imaged, control seen).
pass_metrics <- function(zones, obs = cell_obs, min_ctl_cov = 0) {
  z <- merge(obs, zones, by = "cell")
  m <- z[, .(n_obs = .N, occ = mean(oil_m2 > 0), oil_frac = sum(oil_m2) / sum(cell_m2)), by = .(pass, zone)]
  m <- merge(m, zones[, .(n_zone = .N), by = zone], by = "zone")[, cov := n_obs / n_zone]
  w <- dcast(m, pass ~ zone, value.var = c("occ", "oil_frac", "cov"))
  w <- w[cov_Impact >= min_imp_cov & !is.na(cov_Control) & cov_Control >= min_ctl_cov]
  merge(w, passes, by = "pass")
}

# Plain BACI: before/after means in each zone, difference-in-differences, ratio of ratios.
# y is "occ" or "oil_frac".
did <- function(w, y) {
  I <- w[[paste0(y, "_Impact")]]; C <- w[[paste0(y, "_Control")]]; before <- w$t < restart
  I_b <- mean(I[before]); I_a <- mean(I[!before]); C_b <- mean(C[before]); C_a <- mean(C[!before])
  c(I_before = I_b, I_after = I_a, C_before = C_b, C_after = C_a,
    DiD = (I_a - I_b) - (C_a - C_b), ratio_of_ratios = (I_a / I_b) / (C_a / C_b))
}

# Season-adjusted BACI: regress (impact - control) per pass on an After indicator
# plus calendar-month dummies. Returns the After coefficient.
adj_did <- function(w, y, post_date = restart) {
  gap  <- w[[paste0(y, "_Impact")]] - w[[paste0(y, "_Control")]]
  post <- as.numeric(w$t >= post_date)
  moy  <- factor(month(w$t))
  if (length(unique(post)) < 2) return(NA_real_)
  unname(coef(lm(gap ~ post + moy))["post"])
}

# Block bootstrap: resample whole calendar months (separately for Before and After)
# and recompute `statistic`. Passes in the same month share weather, so they move together.
block_bootstrap <- function(w, statistic, B = n_boot, seed = 1) {
  set.seed(seed)
  months <- split(w, w$month)
  before <- sapply(months, function(m) m$t[1] < restart)
  replicate(B, statistic(rbindlist(months[c(sample(which(before), replace = TRUE),
                                            sample(which(!before), replace = TRUE))])))
}
q_lo <- function(x) unname(quantile(x, .025, na.rm = TRUE))
q_hi <- function(x) unname(quantile(x, .975, na.rm = TRUE))

4 Results

4.1 Where slicks occur

Figure 2: Share of passes with oil in each 1 km cell, before vs after, and the change. The dashed ring is the impact zone around Heritage, Harmony and Hondo (west to east).

The map shows the Channel’s known natural-seep hotspots, Coal Oil Point (~119.9°W) and the Summerland/Carpinteria area (~119.6°W), plus a hotspot near Anacapa and Oxnard in the southeast. In the After period, oil increases along the entire mainland nearshore band, from Hondo east to Carpinteria, not only around the platforms. Surface currents along the mainland coast usually run westward, so seep slicks from Coal Oil Point can drift into the Hondo/Harmony zone.

4.2 Time series

# Main BACI comparison: platform zone (<= R_imp km) vs rest of Channel (> R_excl km)
zones_main <- zone_map(zone_of(grid$dist_km))
per_pass   <- pass_metrics(zones_main)

main <- map_dfr(c("occ", "oil_frac"), function(y) {
  est  <- did(per_pass, y)
  boot <- block_bootstrap(per_pass, function(x) did(x, y)[c("DiD", "ratio_of_ratios")], seed = 1)
  adj  <- adj_did(per_pass, y)
  boot_adj <- block_bootstrap(per_pass, function(x) adj_did(x, y), seed = 2)
  tibble(metric = y, n_passes = nrow(per_pass),
         I_before = est[["I_before"]], I_after = est[["I_after"]], C_before = est[["C_before"]], C_after = est[["C_after"]],
         DiD = est[["DiD"]], DiD_lo = q_lo(boot[1, ]), DiD_hi = q_hi(boot[1, ]),
         rr = est[["ratio_of_ratios"]], rr_lo = q_lo(boot[2, ]), rr_hi = q_hi(boot[2, ]),
         adj_DiD = adj, adj_lo = q_lo(boot_adj), adj_hi = q_hi(boot_adj),
         adj_pct_of_baseline = 100 * adj / est[["I_before"]])
})

# Quarterly means, for the time-series figure
quarterly <- per_pass[, .(n = .N, I = mean(occ_Impact), C = mean(occ_Control)), by = .(q = floor_date(t, "quarter"))][order(q)]

Figure 3: Quarterly slick occurrence in the platform zone vs the rest of the Channel (top), and their difference (bottom).

Before the restart, the platform zone tracked the rest of the Channel, with no trend in the difference. After the restart, it was well above the Channel in Q3–Q4 2025 and again in Q3 2026. In between, in the first half of 2026, it dropped back to baseline, which includes the months right after the full restart. That dip doesn’t fit a simple “more production → more slicks” story. Summer and fall peaks also appear in the pre-period (2023), which is why the seasonal adjustment matters.

4.3 Main BACI estimates

Table 2: Mean per-pass rates (% of imaged area). DiD in percentage points; bootstrap 95% CI in brackets.

Outcome	Impact, before	Impact, after	Control, before	Control, after	DiD	Ratio of ratios	Season-adjusted DiD	Placebo p
Occurrence	0.557%	1.530%	0.720%	0.891%	+0.802 [0.026 to 1.683]	2.22 [1.03–5.46]	+0.661 [-0.069 to 1.430]	0.08
Oiled fraction	0.070%	0.346%	0.166%	0.213%	+0.229 [0.017 to 0.475]	3.84 [1.42–11.00]	+0.204 [-0.003 to 0.456]	0.04

The platform zone went up by a factor of 2.7 and the control by 1.2. That gives a ratio of ratios of 2.22 for occurrence. The oiled-fraction measure shows a larger relative change, because slicks touching the zone after the restart were also bigger (see Slick-level checks). Once the seasonal mismatch is removed, both DiD intervals include zero.

4.4 How unusual is the change?

# Placebo zones: shift the whole platform layout around the Channel on a 4 km lattice
# and re-estimate the season-adjusted effect at each position. If many placebo
# positions show an effect as big as the real one, the real one isn't unusual.
pl_xy   <- st_coordinates(platforms)
cell_xy <- cbind(grid$x, grid$y)
n_imp_real <- zones_main[zone == "Impact", .N]
gb <- st_bbox(grid); lattice <- function(a, b) seq(ceiling(a / 4e3) * 4e3, b, 4e3)
offsets <- expand.grid(dx = lattice(gb[["xmin"]] - min(pl_xy[, 1]), gb[["xmax"]] - max(pl_xy[, 1])),
                       dy = lattice(gb[["ymin"]] - min(pl_xy[, 2]), gb[["ymax"]] - max(pl_xy[, 2])))

placebo <- map_dfr(seq_len(nrow(offsets)), function(i) {
  shifted <- sweep(pl_xy, 2, c(offsets$dx[i], offsets$dy[i]), "+")
  dmin <- apply(sapply(seq_len(nrow(shifted)), function(k)
            sqrt((cell_xy[, 1] - shifted[k, 1])^2 + (cell_xy[, 2] - shifted[k, 2])^2)), 1, min) / 1000
  imp <- dmin <= R_imp
  if (sum(imp) < 0.9 * n_imp_real) return(NULL)       # mostly over land or outside the box
  if (any(grid$dist_km[imp] <= R_excl)) return(NULL)  # overlaps the real platforms
  zones <- zone_map(fifelse(imp, "Impact", fifelse(dmin > R_excl & grid$dist_km > R_excl, "Control", NA_character_)))
  w <- pass_metrics(zones)
  if (nrow(w) < 100) return(NULL)
  tibble(dx = offsets$dx[i], dy = offsets$dy[i], n_passes = nrow(w),
         occ = adj_did(w, "occ"), oil_frac = adj_did(w, "oil_frac"))
})
# Placebo p-value: share of placebo zones with an effect at least as large as the real one
main$placebo_p <- c(mean(placebo$occ >= main$adj_DiD[1]), mean(placebo$oil_frac >= main$adj_DiD[2]))

# Placebo restart date: pretend the restart happened a year earlier, using only pre-restart data
fake_restart <- restart - 365 * 86400
placebo_time <- sapply(c(occ = "occ", oil_frac = "oil_frac"), function(y) adj_did(per_pass[t < restart], y, fake_restart))

# Handy one-row summaries used in the text
occ  <- main %>% filter(metric == "occ")
frac <- main %>% filter(metric == "oil_frac")

Figure 4: Distribution of the season-adjusted BACI effect when the three-platform layout is moved to other ocean locations. The orange line is the actual platform zone.

Of 260 placebo zones, 8% show an occurrence effect at least as large as the real one (4% for oiled fraction). The placebos with effects this large are all over known hotspots: the Coal Oil Point–Carpinteria nearshore seep band, or the hotspot east of Anacapa Island. With the fake restart one year early, the estimated effect is +0.12 points, about 18% of the real estimate. Some of the apparent effect is therefore within the normal year-to-year noise.

4.5 Choice of control area

# Re-run the season-adjusted BACI against alternative control areas. The platforms sit
# 9-13 km off the mainland, so a nearshore control may be a fairer comparison.
not_near <- grid$dist_km > R_excl
nearshore <- not_near & grid$coast_km <= nearshore_km
controls <- setNames(
  list(not_near, nearshore, nearshore & grid$lon < min(platforms_ll$lon), nearshore & grid$lon > max(platforms_ll$lon)),
  c(sprintf("Whole Channel (>%g km)", R_excl), sprintf("Mainland nearshore band (<=%g km of coast)", nearshore_km),
    "Nearshore, west of platforms", "Nearshore, east (Coal Oil Pt-Carpinteria)"))

control_choice <- imap_dfr(controls, function(is_ctl, name) {
  zones <- zone_map(fifelse(grid$dist_km <= R_imp, "Impact", fifelse(is_ctl, "Control", NA_character_)))
  w <- pass_metrics(zones, min_ctl_cov = 0.5)   # control area must be at least half imaged
  map_dfr(c("occ", "oil_frac"), function(y) {
    e <- did(w, y); b <- block_bootstrap(w, function(x) adj_did(x, y), seed = 2)
    tibble(control = name, metric = y, n_passes = nrow(w), n_ctl_cells = sum(is_ctl),
           I_before = e[["I_before"]], I_after = e[["I_after"]], C_before = e[["C_before"]], C_after = e[["C_after"]],
           adj_DiD = adj_did(w, y), lo = q_lo(b), hi = q_hi(b))
  })
})

Table 3: Season-adjusted DiD in occurrence (percentage points) against alternative control areas. ‘Nearshore’ = close to the mainland coast (see settings).

Control	Control cells	Impact before→after	Control before→after	Adjusted DiD [95% CI]
Whole Channel (>10 km)	7833	0.56% → 1.53%	0.72% → 0.89%	+0.66 [-0.07 to 1.43]
Mainland nearshore band (<=12 km of coast)	1671	0.56% → 1.53%	1.29% → 1.98%	+0.22 [-0.63 to 1.01]
Nearshore, west of platforms	328	0.56% → 1.53%	0.26% → 0.23%	+0.76 [0.06 to 1.49]
Nearshore, east (Coal Oil Pt-Carpinteria)	1321	0.56% → 1.53%	1.56% → 2.42%	+0.09 [-0.83 to 0.93]

The platforms are 9–13 km offshore of the mainland, so the whole-Channel control, which includes open water and the island coasts, may not be a fair comparison. The results depend on which control area is used:

• Nearshore west of the platforms (toward Point Conception, few seeps): there’s no change there, so the effect is largest and its interval excludes zero.
• Nearshore east (Coal Oil Point to Carpinteria seeps): this area grew by about the same absolute amount as the platform zone, so the effect is near zero.

In relative terms, though, the platform zone grew more than the seep band did (2.7× vs 1.6×). So one question matters a lot here: would an unaffected platform zone have gained the same number of oiled cells as the seep band, or grown by the same percentage?

4.6 By platform

# Split the impact zone by nearest platform and compare each part with the control zone
per_platform <- map_dfr(sort(platforms$name), function(p) {
  zones <- zone_map(fifelse(grid$dist_km <= R_imp & grid$nearest == p, "Impact",
                            fifelse(grid$dist_km > R_excl, "Control", NA_character_)))
  pass_metrics(zones)[, .(platform = p, n_passes = .N, occ_I = mean(occ_Impact), occ_C = mean(occ_Control)), by = period]
}) %>% mutate(period = factor(period, period_levels)) %>% arrange(platform, period)

# Share of passes with ANY oil near each platform, and whether a single pass dominates
passes_with_oil <- merge(cell_obs, data.table(cell = grid$cell, dist = grid$dist_km, platform = grid$nearest), by = "cell")[dist <= R_imp] %>%
  merge(per_pass[, .(pass, period)], by = "pass") %>%
  .[, .(oiled = any(oil_m2 > 0), cells = sum(oil_m2 > 0)), by = .(pass, period, platform)] %>%
  .[, .(passes = .N, passes_with_oil = sum(oiled), pct_passes_with_oil = round(100 * mean(oiled), 1),
        top_pass_share = round(max(cells) / max(1, sum(cells)), 2)), by = .(platform, period)] %>%
  as_tibble() %>% mutate(period = factor(period, period_levels)) %>% arrange(platform, period)

Figure 5: Occurrence within 5 km of each platform (cells assigned to the nearest platform), by period.

Table 4: Share of passes with any slick in each platform’s zone, and the largest single pass’s share of all oiled cells (a check on whether a few big events drive the result).

Platform	Period	Passes	Passes with oil	% of passes	Largest pass share
Harmony	Before	141	6	4.3%	0.39
	After: phase 1	103	15	14.6%	0.25
	After: phase 2	68	7	10.3%	0.53
Heritage	Before	141	7	5.0%	0.33
	After: phase 1	103	4	3.9%	0.63
	After: phase 2	68	3	4.4%	0.52
Hondo	Before	141	11	7.8%	0.22
	After: phase 1	103	17	16.5%	0.19
	After: phase 2	68	13	19.1%	0.28

• Harmony, the platform that actually restarted in May 2025, has the clearest signal. The share of passes with oil nearby went from 4.3% to 14.6% in phase 1, spread over many passes.
• Hondo also increased (7.8% → 16.5%), even though it wasn’t expected to produce until mid-2026. Hondo and Harmony are only ~4.6 km apart, so their zones touch. Hondo is also the platform nearest the Coal Oil Point seeps, which lie upstream along the coast. Either spillover from Harmony or seep drift could explain this.
• Heritage shows no change, including in phase 2 after it restarted.

4.7 Radius sensitivity

# Does the answer depend on the impact-zone radius? Control = beyond max(R_excl, 2 x radius).
radius <- map_dfr(c(2, 3, 5, 8, 10), function(r) {
  w <- pass_metrics(zone_map(zone_of(grid$dist_km, r, max(R_excl, 2 * r))))
  map_dfr(c("occ", "oil_frac"), function(y) {
    e <- did(w, y)
    b <- block_bootstrap(w, function(x) did(x, y)[["ratio_of_ratios"]], B = 1000, seed = 1)  # fewer draws: 10 runs
    tibble(radius_km = r, metric = y, n_passes = nrow(w), I_before = e[["I_before"]], I_after = e[["I_after"]],
           rr = e[["ratio_of_ratios"]], rr_lo = q_lo(b), rr_hi = q_hi(b))
  })
})

Table 5: Ratio of ratios (unadjusted) for different impact-zone radii. Control = beyond max(10 km, 2 × radius).

Radius (km)	Outcome	Impact before→after	Ratio of ratios [95% CI]
2	Occurrence	0.863% → 2.102%	1.97 [0.78–6.12]
2	Oiled fraction	0.111% → 0.503%	3.51 [0.98–21.38]
3	Occurrence	0.674% → 1.820%	2.18 [0.90–6.28]
3	Oiled fraction	0.079% → 0.437%	4.32 [1.32–17.55]
5	Occurrence	0.557% → 1.530%	2.22 [1.06–5.36]
5	Oiled fraction	0.070% → 0.346%	3.84 [1.46–10.70]
8	Occurrence	0.597% → 1.472%	2.03 [1.10–3.90]
8	Oiled fraction	0.093% → 0.350%	3.00 [1.37–7.45]
10	Occurrence	0.631% → 1.466%	1.92 [1.15–3.27]
10	Oiled fraction	0.107% → 0.358%	2.72 [1.35–5.76]

The relative increase is about 2× for occurrence at every radius from 2 to 10 km. Intervals get narrower as the zone grows, because more cells are averaged, but the bigger zones also take in more of the nearby seep band.

4.8 Slick-level checks

# ---- Individual slicks touching the impact zone, on passes used in the BACI
imp_cells  <- zones_main[zone == "Impact", cell]
imp_slicks <- slick_dt[id %in% slick_cells[cell %in% imp_cells, id] & pass %in% per_pass$pass]
slick_counts <- merge(per_pass[, .(n_passes = .N), by = period],
                      imp_slicks[, .(n_slicks = .N, med_area_km2 = median(area) / 1e6, mean_area_km2 = mean(area) / 1e6), by = period],
                      by = "period", all.x = TRUE)[, per_100 := 100 * n_slicks / n_passes]
area_test <- wilcox.test(area ~ (period == "Before"), data = imp_slicks)   # were After slicks bigger?

# ---- Slicks Cerulean's source model links to each platform (infra_id at any rank)
attribution <- map_dfr(seq_len(nrow(platforms_ll)), function(i) {
  hits <- slick_dt[str_detect(coalesce(source_type_2_ids, ""), paste0("'", platforms_ll$infra_id[i], "'"))]
  merge(passes[, .(n_passes = .N), by = period], hits[, .(n_attr = .N), by = period], by = "period", all.x = TRUE)[
    , platform := platforms_ll$name[i]]
}) %>% as_tibble() %>% mutate(n_attr = replace_na(n_attr, 0L), per_100 = 100 * n_attr / n_passes)

# ---- Main estimate again, ignoring slicks Cerulean classifies as NATURAL seeps (cls 3)
oil_not_natural <- merge(slick_cells[id %in% slick_dt[cls != 3, id]], scene_dt[, .(scene, pass)], by = "scene")[
  , .(oil = sum(oil_m2)), by = .(pass, scene, cell)][, .(oil = max(oil)), by = .(pass, cell)]
obs_not_natural <- merge(cell_obs[, .(pass, cell, cell_m2)], oil_not_natural, by = c("pass", "cell"), all.x = TRUE)[
  , oil_m2 := fcoalesce(oil, 0)]
per_pass_nn <- pass_metrics(zones_main, obs = obs_not_natural)
not_natural <- map_dfr(c("occ", "oil_frac"), function(y) {
  b <- block_bootstrap(per_pass_nn, function(x) adj_did(x, y), B = 1000, seed = 2)
  tibble(metric = y, adj_DiD = adj_did(per_pass_nn, y), lo = q_lo(b), hi = q_hi(b))
})

Table 6: Distinct slicks touching the impact zone, on passes that fully image it.

Period	Passes	Slicks	Per 100 passes	Median area (km²)	Mean area (km²)
Before	141	22	15.6	3.4	7.1
After: phase 1	103	24	23.3	10.4	17.8
After: phase 2	68	19	27.9	4.4	9.7

Slicks touching the platform zone after the restart were larger than before (Wilcoxon rank-sum p = 0.019). Of those After-period slicks, 10 of 43 are labeled natural by Cerulean, compared with 2 of 22 before. Dropping natural-labeled slicks shrinks the season-adjusted occurrence effect to +0.46 points (95% CI -0.11 to 1.15).

Table 7: Slicks Cerulean’s source-attribution model links to each platform’s infrastructure ID (any rank), per 100 passes over the study area.

Platform	Period	Attributed slicks	Per 100 passes
Harmony	Before	2	0.6
	After: phase 1	5	1.9
	After: phase 2	2	1.1
Heritage	Before	6	1.7
	After: phase 1	2	0.8
	After: phase 2	3	1.7
Hondo	Before	6	1.7
	After: phase 1	12	4.7
	After: phase 2	8	4.5

Cerulean’s own attribution agrees in direction. Slicks linked to Hondo and Harmony increased after May 2025, and those linked to Heritage did not. The model attributes a slick to infrastructure based on distance and shape, so it can’t separate a platform discharge from a seep slick that drifts past a platform.

5 Interpretation

• There is an increase near the platforms that shouldn’t be dismissed. It lines up in time and place with the one platform that restarted first (Harmony), it’s spread across many passes rather than one or two big events, and Cerulean’s attribution model points the same way.
• The evidence is not strong enough to attribute the increase to the restart.
  1. After adjusting for season, the effect’s confidence interval includes zero, and about 8% of random placebo zones show a change this large.
  2. The nearby seep band got oilier by a similar absolute amount over the same period, which suggests something affected the whole nearshore strip: wind, currents, seep activity, or Cerulean model changes.
  3. The timing is inconsistent with production. The platform zone returned to baseline in early 2026, right after the full restart.
  4. Heritage, which restarted in April 2026, shows nothing, and Hondo rose before its own restart.
• Bottom line: the data are consistent with a modest increase in slicks around Harmony/Hondo after May 2025. They are about equally consistent with a nearshore-wide increase in seep slicks drifting past the platforms. More post-restart data, especially after Hondo comes online, should help separate the two.

6 Limitations

• Detectability. SAR only shows slicks in moderate wind (roughly 2–10 m/s). Wind isn’t controlled for directly; comparing zones on the same pass removes Channel-wide wind conditions, but not local differences.
• Model changes. Cerulean’s detection model has changed over time. That’s only absorbed if the changes affected the whole Channel equally.
• Seeps vs. platforms. Seep slicks can’t be reliably separated from platform slicks. Cerulean’s class labels are model output, not ground truth; only 8 slicks have human review.
• Coarse coastline. Coastline is Natural Earth 1:50m, which is coarse near shore and around the islands.
• One treated site. There’s a single treated area, so conventional standard errors aren’t meaningful. The placebo tests are the more honest measure of uncertainty.
• Short After period (May 2025–Sep 2026), which runs two summers against one winter.

7 Next steps

1. Add wind speed (ERA5 or NDBC buoys 46053/46054) as a pass-level covariate, or restrict to passes with good detection conditions.
2. Use slick shape and orientation: slicks from a platform should start at the platform and stretch downwind or downcurrent, while seep slicks start over known seep fields.
3. Re-run as more post-restart data accumulates, especially after Hondo comes online (see Updating this report).
4. Compare against production volumes by platform, if monthly data are available (BSEE OGOR).

8 Updating this report

Everything (download, analysis, figures and this text) comes from one file, report.qmd. Open Show code under any section to see how that step was done, or use the </> Code menu at the top to show all code or view the full source.

• Change a setting (dates, restart dates, platforms, zone sizes): edit the settings block at the top of report.qmd. The settings used for this version are listed below.
• Get the latest satellite data: set refresh_data <- TRUE, render once, set it back to FALSE, and commit the updated data/ folder.
• Render: quarto render report.qmd (about 2 minutes), or click Render in RStudio. On GitHub, any change pushed to report.qmd or data/ re-renders and republishes the report automatically.

Numbers in the text update automatically. The written interpretation was written for data through September 2026, so reread it after big changes.

8.1 Settings used

Data run from Jan 1, 2023 to Sep 26, 2026.

# ============================================================================
# SETTINGS: everything you might want to change lives here.
# ============================================================================

# ---- Study area -------------------------------------------------------------
# Bounding box in longitude/latitude (WGS84). Used both for the data download
# and as the analysis area (the ocean inside this box).
bbox <- c(xmin = -120.65, ymin = 33.85, xmax = -119.0, ymax = 34.50)

# ---- Dates ------------------------------------------------------------------
# Period of slicks and satellite scenes to use (inclusive, UTC). An end date in
# the future is fine: you get everything available up to today.
date_start <- "2023-01-01"
date_end   <- "2026-12-31"

# ---- Restart dates ----------------------------------------------------------
# "Before" = everything before `restart`. "After" is split into two phases.
restart  <- as.POSIXct("2025-05-01", tz = "UTC")  # Harmony test restart: start of "After"
restart2 <- as.POSIXct("2026-03-14", tz = "UTC")  # shipping resumes; start of "After: phase 2"

# ---- Platforms --------------------------------------------------------------
# Locations from the BOEM Santa Ynez Unit Development and Production Plan.
# infra_id = the platform's ID in Cerulean's infrastructure database, used to
# count slicks Cerulean's model attributes to each platform.
platforms_ll <- data.frame(
  name     = c("Hondo",     "Harmony",   "Heritage"),
  lon      = c(-120.1203,   -120.1675,   -120.2786),
  lat      = c(34.3907,     34.3767,     34.3502),
  infra_id = c("74569",     "117006",    "430365")
)

# ---- Zones ------------------------------------------------------------------
R_imp        <- 5     # km. Impact zone = ocean within this distance of any platform
R_excl       <- 10    # km. Control zone = ocean farther than this from every platform
cell_m       <- 1000  # m. Size of the square grid cells the ocean is divided into
min_imp_cov  <- 0.9   # a satellite pass is used only if it images at least this share of the impact zone
nearshore_km <- 12    # km. "Nearshore" control areas = within this distance of the mainland coast

# ---- Statistics -------------------------------------------------------------
n_boot <- 2000        # bootstrap draws for 95% intervals (more = slower but more stable)

# ---- Data -------------------------------------------------------------------
# FALSE: use the snapshot saved in data/ (fast, reproducible).
# TRUE:  re-download everything from the Cerulean API (a few minutes; no key needed).
refresh_data <- FALSE
cerulean_api <- "https://api.cerulean.skytruth.org/collections"

9 Sources

• SkyTruth Cerulean API: https://api.cerulean.skytruth.org
• World Oil (May 19, 2025): Sable Offshore restarts oil production from Santa Ynez Unit
• Noozhawk: Sable restarts oil production, transport after federal order
• Oil & Gas Journal (Apr 2, 2026): Sable set to resume production at second platform
• BOEM: Santa Ynez Unit Development and Production Plan (platform locations)