Local DriftCaster: trajectory calculations, ensembles, and updates

This document describes the equations and operational workflow implemented by the Local DriftCaster trajectory viewer. It is intended to make every displayed track reproducible from its archived forcing inputs and solver settings.

The active viewer contains two kinds of cases:

nine observation-backed validation cases for Hooper Bay, Englullrarmiut, and Freebird; and
five hypothetical Okinawa priority forecasts, each with a five-point north-south release line and no observations yet.

The operational calculation code is in ../data/SPOTTER/driftcast-local-forcing-qc-2026-08-18/. The generated viewer and its currently published case files are in spotter_validation_viewer/.

1. Coordinate and notation conventions

All horizontal vectors use east-north order:

[ \mathbf U=(u,v), ]

where (u) is eastward and (v) is northward. Velocities are in (\mathrm{m,s^{-1}}), positions are longitude/latitude in decimal degrees, and internal times are UTC Unix seconds.

Symbol	Meaning	Source variable(s)
(\mathbf U_c)	Eulerian surface-ocean current	SMOC `uo`, `vo` at the shallowest depth
(\mathbf U_S)	surface Stokes drift	SMOC `vsdx`, `vsdy` for Okinawa; wave-product `VSDX`, `VSDY` for validation
(\overline{\mathbf U}_S)	draft-averaged Stokes drift used by the validation cases	derived from `VSDX`, `VSDY`, and `VTPK`
(\mathbf U_{10})	10 m wind	IFS `u10`, `v10`
(\mathbf U_T)	tidal drift	SMOC `utide`, `vtide`
(\mathbf U_{\rm total})	full SMOC surface velocity	SMOC `utotal`, `vtotal`
(\mathbf V)	modeled buoy velocity	calculated
(\mathbf x=(x,y))	local east/north position	calculated
((\lambda,\phi))	longitude and latitude	calculated or observed
(R)	mean Earth radius	6,371,000 m

The full SMOC variables are treated as the composite velocity supplied by the product. The code does not add Stokes drift or tidal drift to (\mathbf U_{\rm total}) again.

2. Forcing products

Operational Okinawa forecasts

The Okinawa calculations use:

Copernicus Marine SMOC dataset cmems_mod_glo_phy_anfc_merged-uv_PT1H-i, including uo, vo, vsdx, vsdy, utide, vtide, utotal, and vtotal;
ECMWF deterministic IFS 10 m wind, u10 and v10; and
the Copernicus global wave forecast for viewer-only wave height, period, and direction backgrounds.

All velocity components used by a single Okinawa run come from one archived SMOC snapshot. This prevents the simple and full-SMOC tracks from silently mixing forecast cycles.

The regional SMOC request spans 115-145 degrees E and 15-35 degrees N, from one hour before release through one hour after the 72-hour forecast. The extra hour on each side is required by the centered material-derivative calculation.

Observation-backed validation cases

The nine fixed examples use the locally archived currents, waves, and winds associated with their historical release windows. Wind switches from ERA5 to the archived IFS forecast. If a gap exists between the last ERA5 time and first IFS time, the code linearly bridges the two endpoints.

The validation cases use a draft-averaged Stokes velocity. With peak period (T_p), deep-water wave number

[ k=\frac{(2\pi/T_p)^2}{g}, ]

draft (d=0.155\ \mathrm{m}), and (z=2kd), the factor is

[ F_d=\frac{1-e^{-z}}{z}. ]

Therefore

[ \overline{\mathbf U}_S=F_d\mathbf U_S. ]

The Okinawa forecast generator instead uses the SMOC vsdx/vsdy surface values directly. This is an intentional implementation difference that must be retained when reproducing the archived results.

3. Interpolation

Forcing is sampled at every ensemble member’s current time and position. The implementation performs:

bilinear interpolation over the four surrounding longitude/latitude cells;
renormalization of the spatial weights if one or more corners are missing;
linear interpolation between the two surrounding forcing times.

For a spatial field value (q), the four bilinear weights are

[ \begin{aligned} w*{00}&=(1-f_y)(1-f_x), & w*{01}&=(1-fy)f_x,\ w{10}&=fy(1-f_x), & w{11}&=f_yf_x, \end{aligned} ]

and the missing-value-aware estimate is

[ qs=\frac{\sum{ij\in\mathcal F}w*{ij}q*{ij}} {\sum*{ij\in\mathcal F}w*{ij}}, ]

where (\mathcal F) contains the finite corners. Linear time interpolation is

[ q(t)=(1-f_t)q_s(t_0)+f_tq_s(t_1). ]

Coordinates outside an input axis are clamped to the nearest boundary by the interpolator. Operational input coverage is therefore checked before a run so time-boundary clamping is not used as a substitute for missing forecast data.

4. Trajectory models

4.1 Standard simple formulation

The general weighted model is

[ \boxed{\mathbf V=\alphac\mathbf U_c+ \alpha_s\mathbf U_S+ \alpha_w\mathbf U{10}}. ]

The active standard coefficients are

Coefficient	Value	Interpretation
(\alpha_c)	1.0	full current
(\alpha_s)	1.0	full Stokes term
(\alpha_w)	0.012	1.2% of 10 m wind

Thus the active model is

[ \boxed{\mathbf V0=\mathbf U_c+\mathbf U_S+0.012\mathbf U{10}}. ]

The windage coefficient is (\alpha_w=0.012), or 1.2% of the 10 m wind.

For the validation examples, (\mathbf U_S) in this expression is replaced by the draft-averaged (\overline{\mathbf U}_S) defined above.

The simple and full-SMOC models use a 60 s integration step. For a velocity held over one step (\Delta t), the spherical position update is

[ \phi*{n+1}=\phi_n+\frac{v_n\Delta t}{R}, \qquad \lambda*{n+1}=\lambdan+ \frac{u_n\Delta t}{R\cos\phi{n+1}}, ]

with angular increments converted from radians to degrees. Positions are saved hourly, including the release position, producing 73 points for a 72-hour run.

4.2 Full SMOC plus windage

The third model shown for the Okinawa cases and the existing validation cases is

[ \boxed{\mathbf V*{\rm SMOC}=\mathbf U*{\rm total}+ 0.012\mathbf U_{10}}. ]

utotal and vtotal are used exactly once. Although tidal and Stokes fields are available as viewer backgrounds, they are not explicitly added to this trajectory because they are already represented by the full SMOC product.

4.3 Finite-inertia formulation

The finite-inertia model evolves both position and buoy velocity. Define the water velocity and relative velocities as

[ \mathbf UL=\mathbf U_c+\mathbf U_S, \qquad \mathbf q_w=\mathbf U_L-\mathbf V, \qquad \mathbf q_a=\mathbf U{10}-\mathbf V. ]

Again, the validation cases use (\overline{\mathbf U}_S), while Okinawa uses the SMOC surface Stokes velocity.

The effective mass includes added mass:

[ M=mb+C_A\rho_wV{sub}. ]

For each fluid (j\in{w,a}),

[ Re_j=\frac{L_j|\mathbf q_j|}{\nu_j}, ]

and the drag coefficient is evaluated in three Reynolds-number regimes:

[ C_D(Re)=0 \qquad \text{for } Re\le10^{-12}, ]

[ C_D(Re)=\frac{24}{Re}\left(1+0.15Re^{0.687}\right) \qquad \text{for } 0<Re<1000, ]

[ C_D(Re)=0.44 \qquad \text{for } Re\ge1000. ]

The water and air drag accelerations are represented by

[ k_w=\frac{\tfrac12\rho_w C_D(Re_w)A_w|\mathbf q_w|}{M}, \qquad k_a=\frac{\tfrac12\rho_a C_D(Re_a)A_a|\mathbf q_a|}{M}. ]

With Coriolis parameter

[ f=2\Omega\sin\phi, ]

the exact component equations implemented in the solver are

[ \boxed{ \begin{aligned} \dot Vx &= \frac{D U{L,x}}{Dt} +kwq{w,x}+kaq{a,x}-fq*{w,y},\ \dot V_y &= \frac{D U*{L,y}}{Dt} +kwq{w,y}+kaq{a,y}+fq_{w,x}. \end{aligned}} ]

The water material derivative is

[ \frac{D\mathbf UL}{Dt}= \frac{\partial\mathbf U_L}{\partial t} +U{L,x}\frac{\partial\mathbf UL}{\partial x} +U{L,y}\frac{\partial\mathbf U_L}{\partial y}. ]

It is evaluated with centered differences using (\delta t=1800\ \mathrm{s}) and (\delta x=\delta y=10{,}000\ \mathrm{m}):

[ \frac{\partial\mathbf U_L}{\partial t}\approx \frac{\mathbf U_L(t+\delta t)-\mathbf U_L(t-\delta t)}{2\delta t}, ]

with analogous centered differences in east and north. The derivative is refreshed every 600 s and the local water and wind velocities every 60 s.

The solver uses forward Euler:

[ \mathbf x*{n+1}=\mathbf x_n+\mathbf V_n\Delta t, \qquad \mathbf V*{n+1}=\mathbf V_n+\dot{\mathbf V}_n\Delta t. ]

The Okinawa finite-inertia runs use (\Delta t=1\ \mathrm{s}). The northern validation cases use 4 s and the Freebird validation cases use 1 s. Okinawa is initialized with the local standard simple-model velocity. The validation cases instead initialize velocity from a linear fit to the observed track over approximately +/-30 minutes for the northern buoys and +/-15 minutes for Freebird.

Finite-inertia constants

Variable	Value	Units	Meaning
(m_b)	7.45	kg	buoy mass
(V_{sub})	0.0072683	m3	submerged volume
(C_A)	0.5	dimensionless	added-mass coefficient
(M)	11.17500375	kg	effective mass
(A_w)	0.03643	m2	submerged frontal area
(L_w)	0.2993	m	water-side characteristic length
(A_a)	0.0406875	m2	exposed frontal area
(L_a)	0.42	m	air-side characteristic length
(\rho_w)	1025	kg m-3	seawater density
(\rho_a)	1.225	kg m-3	air density
(\nu_w)	1.05e-6	m2 s-1	water kinematic viscosity
(\nu_a)	1.50e-5	m2 s-1	air kinematic viscosity
(\Omega)	7.2921159e-5	s-1	Earth rotation rate
(g)	9.80665	m s-2	gravitational acceleration
(d)	0.155	m	draft used for validation Stokes averaging
(R)	6,371,000	m	Earth radius

These are fixed engineering parameters, not fitted separately for each case.

5. Ensemble construction

Only the initial position is perturbed. For each member,

[ (\delta x,\delta y)\sim\mathcal N \left(\mathbf0,\sigma^2\mathbf I\right), \qquad \sigma=1000\ \mathrm{m}. ]

The conversion to geographic coordinates is

[ \delta\phi=\frac{\delta y}{R}, \qquad \delta\lambda=\frac{\delta x}{R\cos\phi_0}. ]

The 1 km setting is the standard deviation of each east/north component, not a hard radial limit. For an ideal two-dimensional Gaussian, its initial 80% radial quantile is approximately 1.794 km.

Sampling is antithetic: every random offset ((\delta x,\delta y)) is paired with ((-\delta x,-\delta y)). This centers each finite ensemble exactly on the requested release. The same offsets are shared across models, and the five Okinawa sites also share the same offset set, making model/site comparisons paired rather than confounded by different random draws.

Current Okinawa settings are:

100 members per site per model;
five sites per priority area;
500 tracks per model and 1,500 tracks per priority forecast;
three models: standard simple, finite inertia, and full SMOC plus windage;
random seed 20260825.

Across all five priorities this is 7,500 trajectories per initialization time. The four Okinawa initialization batches contain 30,000 trajectories in total. The nine validation cases contain 100 members per model, with case-specific deterministic seeds derived from 20260819.

6. The 80% trajectory cone

At each saved hour, the viewer product calculates component-wise median longitude and latitude. It then computes every member’s great-circle distance from that median and stores the 0.8 quantile as (r_{80}(t)).

Great-circle distance uses the haversine expression

[ d=2R\arcsin\sqrt{ \sin^2\left(\frac{\Delta\phi}{2}\right)+ \cos\phi_1\cos\phi_2 \sin^2\left(\frac{\Delta\lambda}{2}\right)}. ]

The cone is a single merged geometry, rather than many overlapping translucent circles. In a local tangent plane, the code:

offsets the median centerline by +/-(r_{80}(t)) along its local normal;
joins adjacent left/right offsets with quadrilaterals;
unions those pieces with a circular buffer at every hourly center;
uses 16 segments per quadrant for the round buffers.

This creates a clean outline with a rounded end cap whose radius is the final 80% spread. It is a visualization of hourly marginal 80% spread, not a formal simultaneous 80% confidence region for the entire trajectory.

For a five-site forecast, each release group has its own cone. A combined outline is also stored; disconnected pieces fall back to their convex hull for viewer rendering.

7. Okinawa release geometry and timing

Each priority transect is evaluated at four initialization times. Every case runs for 72 hours and saves positions hourly.

Initialization in JST	Initialization in UTC	Forecast end in UTC
2026-08-25 09:00 JST	2026-08-25 00:00 UTC	2026-08-28 00:00 UTC
2026-08-25 21:00 JST	2026-08-25 12:00 UTC	2026-08-28 12:00 UTC
2026-08-26 09:00 JST	2026-08-26 00:00 UTC	2026-08-29 00:00 UTC
2026-08-26 21:00 JST	2026-08-26 12:00 UTC	2026-08-29 12:00 UTC

The same release geometry, initial perturbations, random seed, coefficients, models, and numerical settings are used at every initialization. Only the release time and therefore the sampled time-evolving forcing change. The viewer groups the twenty priority cases by JST initialization time.

Each provided midpoint is site 3. Sites 1-5 are at latitude offsets (+2,+1,0,-1,-2) arcminutes with unchanged longitude, so adjacent nominal release points are one nautical mile apart north-south.

Priority	Center latitude	Center longitude
1	27 deg 16 min 0.61 sec N	127 deg 38 min 41.64 sec E
2	26 deg 55 min 47.44 sec N	127 deg 16 min 57.56 sec E
3	26 deg 12 min 24.34 sec N	126 deg 40 min 54.52 sec E
4	25 deg 59 min 25.80 sec N	126 deg 39 min 44.00 sec E
5	25 deg 59 min 0.06 sec N	127 deg 10 min 45.89 sec E

There are no observations for these hypothetical releases, so no forecast error or skill score is assigned to them before deployment.

8. Validation metrics

For observation-backed cases, pointwise position error is the great-circle distance

[ ex(t)=d\left(\mathbf x{model}(t),\mathbf x_{obs}(t)\right). ]

The stored headline metrics are:

[ E*{mean}=\operatorname{median}*{member} \left[\frac{1}{73}\sum_{h=0}^{72}e_x(h)\right], ]

[ E*{72}=\operatorname{median}{member}\left[e_x(72)\right], \qquad S{72}=r*{80}(72). ]

Hourly interval velocity is reconstructed from consecutive positions:

[ u_h=\frac{R\cos\bar\phi,\Delta\lambda}{3600}, \qquad v_h=\frac{R,\Delta\phi}{3600}. ]

Velocity error is

[ ev=\left|\mathbf V{model}-\mathbf V_{obs}\right|. ]

With observed-motion unit vector (\hat{\mathbf t}=\mathbf V*{obs}/|\mathbf V*{obs}|) and perpendicular (\hat{\mathbf n}=(-\hat t_y,\hat t_x)), the diagnostic magnitudes are

[ e*{along}=|\Delta\mathbf V\cdot\hat{\mathbf t}|, \qquad e*{cross}=|\Delta\mathbf V\cdot\hat{\mathbf n}|. ]

For candidate-minus-reference comparisons, positive improvement is defined as

[ Ix=e{x,0}-e*{x,a}, \qquad I_v=e*{v,0}-e_{v,a}. ]

Thus positive values mean the candidate has less error. Error growth rate is approximated by

[ \dot e(t)\approx\frac{e(t+\Delta t)-e(t)}{\Delta t}. ]

The primary coefficient-tuning objective used in the diagnostic suite is the median across forecasts of mean 0-72 h position error:

[ J(a)=\operatorname{median}{forecast} \left[\frac{1}{T}\int_0^T |\mathbf x_a(t)-\mathbf x{obs}(t)|,dt\right]. ]

9. Coefficient experiments versus active settings

The repository contains earlier tuning products; these do not silently alter the active viewer model.

Experiment	Values tested	Coefficients held fixed
windage	0.000-0.015 in 0.001 increments	(\alpha_c=1,\alpha_s=1)
current	0.5-1.5 in 0.1 increments	(\alpha_s=1,\alpha_w=0.012)
Stokes	0, 0.25, 0.5, …, 2.0	(\alpha_c=0.6,\alpha_w=0.012)
hourly rolling comparison	(\alpha_c=1.0) versus 0.6	(\alpha_s=1,\alpha_w=0.012)

The tuned (\alpha_c=0.6) case remains an analysis experiment. The active ensemble forecast labeled Simple (1, 1, 0.012) uses (\alpha_c=1.0).

10. Forecast selection and automatic updates

The macOS LaunchAgent runs poll_forecasts.py every 900 s (15 minutes) while the user is logged in and the Mac is awake. The workflow is:

poll and archive eligible ECMWF IFS cycles;
poll and archive Copernicus forecast updates;
if both pollers succeed, check whether the Okinawa cases need rebuilding;
archive the complete new run first;
atomically replace the five cases belonging to that initialization and update the manifest without removing the other initialization batches.

IFS cycles are initialized at 00, 06, 12, and 18 UTC. The 00/12 UTC products extend to 240 h; 06/18 UTC extend to 144 h. For the Okinawa experiment the refresh controller selects the latest initialization that is no later than the release time and still covers the complete 72-hour interval:

[ t*{cycle}\le t*{release}, \qquad t*{cycle}+H*{available}\ge t_{release}+72\ \mathrm{h}. ]

This is the closest eligible pre-release forecast, not simply the newest file on disk.

A fresh regional SMOC snapshot is requested no more often than every six hours by default. Its coordinate arrays and all eight trajectory variables are content-hashed with SHA-256. Trajectories are rebuilt only if the eligible IFS cycle changes, the SMOC content hash changes, or a forced run is requested.

Automatic updating stops at the release time. This prevents post-release model updates from introducing hindsight into later validation. Retrospective runs are possible only with the explicit --allow-after-release option and remain timestamped in the archive.

The refresh refuses to run with less than 15 GiB of free disk space, validates full SMOC time coverage and required variables, and uses a filesystem lock to prevent concurrent rebuilds.

Manual update commands

From the repository root:

cd data/SPOTTER/driftcast-local-forcing-qc-2026-08-18
 
# Inspect the selected IFS cycle and whether a SMOC refresh is due.
.venv/bin/python refresh_okinawa_priority_forecasts.py --check-only
 
# Poll all operational sources, then refresh if inputs changed.
.venv/bin/python poll_forecasts.py
 
# Force a new pre-release SMOC snapshot and recalculation.
.venv/bin/python refresh_okinawa_priority_forecasts.py --force

The trajectory generator itself uses the Python interpreter named by DRIFTCAST_TRAJECTORY_PYTHON; if unset, it prefers /Users/wyn/miniforge3/envs/Pacific/bin/python and otherwise falls back to the current interpreter.

11. Archive and provenance

The active viewer contains twenty Okinawa cases: five priorities at four initialization times. Recalculating one initialization replaces only its five active viewer JSON files. Historical calculations are never overwritten in the archive.

Immutable runs live at:

data/SPOTTER/driftcast-local-forcing-qc-2026-08-18/
  forecast_archive/okinawa_priority_forecasts/
    index.json
    runs/
      YYYYMMDDTHHMMSSZ/
        run.json
        cases/
          ensemble_okinawa_priority-1_YYYYMMDDTHHZ.json
          ...
          ensemble_okinawa_priority-5_YYYYMMDDTHHZ.json

Every run records:

UTC calculation time and unique run ID;
release and forecast-end time;
IFS path, initialization cycle, SHA-256, and byte count;
SMOC snapshot path, retrieval time, content SHA-256, byte count, and variable list;
generator path and SHA-256;
ensemble size, position sigma, seed, and finite-inertia time step; and
all five complete case payloads, including every member track, medians, cones, fields, and viewer metadata.

JSON publication is atomic: data are written to a temporary file and moved into place only after serialization succeeds. A failed recalculation therefore does not leave a partially written active case.

For later deployment validation, match observations to archived cases by forecast_run_id, release.time, and release_sites[].id. Do not use only the current viewer JSON, because it intentionally represents the latest eligible run rather than the full forecast history.

To see the calculation date and forcing cycles for the latest archived run:

jq '.runs[-1]' \
  data/SPOTTER/driftcast-local-forcing-qc-2026-08-18/forecast_archive/okinawa_priority_forecasts/index.json

12. Viewer forcing layers

The following time-evolving backgrounds are packaged with each Okinawa case:

Eulerian ocean-current speed and vectors;
Stokes-drift speed and vectors;
significant wave height with mean wave direction;
peak wave period with peak direction;
mean wave period;
10 m wind speed and vectors;
full SMOC total speed and vectors; and
tidal-drift speed and vectors.

Viewer fields are snapshots every 3 h, independent of the solver sampling. The high-resolution immediate context uses a 31 x 21 display grid. A coarser 29 x 15 East China Sea field over 115-145 degrees E and 15-35 degrees N supplies basin-wide context. These display grids do not control trajectory resolution; the trajectories sample the underlying forcing grids directly.

13. Running the viewer