Okinawa Trajectory Calculation Methods

Local DriftCast: trajectory calculations, ensembles, and updates

This document describes the equations and operational workflow implemented by the Local DriftCast 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}{\mathrm{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}}_{\mathrm{total}}$	full SMOC surface velocity	SMOC utotal, vtotal
$\mathbf{V}$	modeled buoy velocity	calculated
$\mathbf{x}=(x,y)$	local east/north position	calculated
$(\lambda ,\varphi )$	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}}_{\mathrm{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 ocean-velocity components used by a single Okinawa run come from one archived SMOC snapshot. This prevents the SMOC and finite-inertia 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:

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

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

$$\begin{array}{rlrl}{w}_{00}& =(1-f_y)(1-f_x),& w_{01}& =(1-f_y)f_x,\\ w_{10}& =f_y(1-f_x),& w_{11}& =f_y{f}_x,\end{array}$$

and the missing-value-aware estimate is

$$q_s=\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_t{q}_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

The active Okinawa viewer compares two trajectory formulations: the SMOC additive model and the finite-inertia model. The earlier standalone simple track is no longer plotted as a separate model because it differs from the SMOC formulation primarily by omitting the tidal contribution.

4.1 SMOC formulation: simple model plus tides

The SMOC trajectory retains the direct additive structure of the simple model and adds tidal drift. In physical-component notation, the displayed model is

$$\boxed{{\mathbf{V}}_{\mathrm{SMOC}}={\mathbf{U}}_c+{\mathbf{U}}_S+{\mathbf{U}}_T+0.012{\mathbf{U}}_{10}}.$$

The active coefficients are

Coefficient	Value	Interpretation
${\alpha }_c$	1.0	full current
${\alpha }_s$	1.0	full Stokes term
${\alpha }_T$	1.0	full tidal term
${\alpha }_w$	0.012	1.2% of 10 m wind

For the operational calculation, SMOC supplies these water-motion components as the composite utotal and vtotal variables. The equation evaluated by the integrator is therefore

$$\boxed{{\mathbf{V}}_{\mathrm{SMOC}}={\mathbf{U}}_{\mathrm{total}}+0.012{\mathbf{U}}_{10}},{\mathbf{U}}_{\mathrm{total}}\equiv {\mathbf{U}}_c+{\mathbf{U}}_S+{\mathbf{U}}_T.$$

This is the former simple current + Stokes + windage formulation with the tidal term included. The windage coefficient remains ${\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.

utotal and vtotal are used exactly once. The separately published current, Stokes, and tidal fields remain available as diagnostic viewer backgrounds; they are not added a second time to the SMOC trajectory.

The SMOC model uses a 60 s integration step. For a velocity held over one step $\Delta t$, the spherical position update is

$${\varphi }_{n+1}={\varphi }_n+\frac{v_n\Delta t}{R},{\lambda }_{n+1}={\lambda }_n+\frac{u_n\Delta t}{R\cos {\varphi }_{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 Finite-inertia formulation

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

$${\mathbf{U}}_L={\mathbf{U}}_c+{\mathbf{U}}_S,{\mathbf{q}}_w={\mathbf{U}}_L-\mathbf{V},{\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=m_b+C_A{\rho }_w{V}_{sub}.$$

For each fluid $j\in \{w,a\}$,

$$R{e}_j=\frac{L_j|{\mathbf{q}}_j|}{{\nu }_j},$$

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

$$C_D(Re)=0\text{for }Re\le 10^{-12},$$ $$C_D(Re)=\frac{24}{Re}\left(1+0.15R{e}^{0.687}\right)\text{for }0<Re<1000,$$ $$C_D(Re)=0.44\text{for }Re\ge 1000.$$

The water and air drag accelerations are represented by

$$k_w=\frac{\frac{1}{2}{\rho }_w{C}_D(R{e}_w)A_w|{\mathbf{q}}_w|}{M},k_a=\frac{\frac{1}{2}{\rho }_a{C}_D(R{e}_a)A_a|{\mathbf{q}}_a|}{M}.$$

With Coriolis parameter

$$f=2\Omega \sin \varphi ,$$

the exact component equations implemented in the solver are

$$\boxed{\begin{array}{rl}{\dot{V}}_x& =\frac{D{U}_{L,x}}{Dt}+k_w{q}_{w,x}+k_a{q}_{a,x}-f{q}_{w,y},\\ {\dot{V}}_y& =\frac{D{U}_{L,y}}{Dt}+k_w{q}_{w,y}+k_a{q}_{a,y}+f{q}_{w,x}.\end{array}}$$

The water material derivative is

$$\frac{D{\mathbf{U}}_L}{Dt}=\frac{\partial {\mathbf{U}}_L}{\partial t}+U_{L,x}\frac{\partial {\mathbf{U}}_L}{\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,{\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 current + Stokes + windage velocity. Its subsequent finite-inertia evolution remains unchanged and does not substitute the SMOC tidal composite for ${\mathbf{U}}_L$. 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(0,{\sigma }^2\mathbf{I}\right),\sigma =1000\mathrm{m}.$$

The conversion to geographic coordinates is

$$\delta \varphi =\frac{\delta y}{R},\delta \lambda =\frac{\delta x}{R\cos {\varphi }_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,000 displayed tracks per priority forecast;
• two displayed models: SMOC (current + Stokes + tides + 0.012 wind) and finite inertia;
• random seed 20260825.

Across all five priorities this is 5,000 displayed trajectories per initialization time. The four Okinawa initialization batches contain 20,000 displayed 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 \varphi }{2}\right)+\cos {\varphi }_1\cos {\varphi }_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:

1. offsets the median centerline by +/-$r_{80}(t)$ along its local normal;
2. joins adjacent left/right offsets with quadrilaterals;
3. unions those pieces with a circular buffer at every hourly center;
4. 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

$$e_x(t)=d\left({\mathbf{x}}_{model}(t),{\mathbf{x}}_{obs}(t)\right).$$

The stored headline metrics are:

$$E_{mean}={\mathrm{median}}_{member}\left[\frac{1}{73}\sum _{h=0}^{72}{e}_x(h)\right],$$ $$E_{72}={\mathrm{median}}_{member}\left[e_x(72)\right],S_{72}=r_{80}(72).$$

Hourly interval velocity is reconstructed from consecutive positions:

$$u_h=\frac{R\cos \bar{\varphi }\Delta \lambda }{3600},v_h=\frac{R\Delta \varphi }{3600}.$$

Velocity error is

$$e_v=\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}}|,e_{cross}=|\Delta \mathbf{V}\cdot \hat{\mathbf{n}}|.$$

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

$$I_x=e_{x,0}-e_{x,a},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)={\mathrm{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 and the earlier no-tide simple track remain analysis products. Neither is presented as a separate active Okinawa model. The active additive track uses unit current, Stokes, and tidal contributions plus $0.012{\mathbf{U}}_{10}$ through the SMOC composite formulation in Section 4.1.

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:

1. poll and archive eligible ECMWF IFS cycles;
2. poll and archive Copernicus forecast updates;
3. if both pollers succeed, check whether the Okinawa cases need rebuilding;
4. archive the complete new run first;
5. 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},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;
• combined SMOC current + Stokes + tidal 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

From the repository root:

./code/start

Open http://localhost:8000/. A local web server is required because browsers do not allow the case files to be fetched reliably from a file:// URL. Keep the terminal process running and press Ctrl-C to stop it.

To use another port:

SPOTTER_VIEWER_PORT=8765 ./code/start

The older observation-validation builder can still be run with:

python3 code/spotter_validation.py

That command regenerates the legacy validation product; it is separate from the archived operational Okinawa refresh workflow described above.

14. Primary implementation files

• ../data/SPOTTER/driftcast-local-forcing-qc-2026-08-18/notebooks/generate_okinawa_line_ensemble_case.py: five-site Okinawa solver and viewer payloads.
• ../data/SPOTTER/driftcast-local-forcing-qc-2026-08-18/notebooks/generate_ensemble_validation_cases.py: nine observation-backed ensembles, finite-inertia equations, and cone geometry.
• ../data/SPOTTER/driftcast-local-forcing-qc-2026-08-18/refresh_okinawa_priority_forecasts.py: forcing selection and recalculation policy.
• ../data/SPOTTER/driftcast-local-forcing-qc-2026-08-18/forecast_run_archive.py: immutable run storage and atomic publishing.
• ../data/SPOTTER/driftcast-local-forcing-qc-2026-08-18/poll_forecasts.py: operational polling entry point.
• spotter_validation_viewer/index.html: Local DriftCast viewer.