simulate_particles() |> simulate_particles(): a 2D stochastic Lagrangian particle simulation using nearest-neighbour advection from gridded tiltmeter flow fields.
Approach: coralseed operates as a settlement and post-settlement modelling framework that can be run downstream of, or alongside, particle dispersal simulations by translating larval movement pathways into probabilistic settlement and demographic outcomes across fine-scale habitat maps. coralseed can also be applied to particle tracks generated from simplified advection assumptions, empirical current-meter observations, or user-defined dispersal scenarios. Field observations, including local flow-meter instrumental records, can be used to parameterise movement by converting measured current speed and direction into lines of movement from release locations over biologically relevant time steps.
Methods: Here we use a simplified, observation-driven 2D Lagrangian particle-tracking model to explore short-term passive dispersal from three release locations. Particles were advected using spatially interpolated logger-derived eastward and northward velocities at 5 min timesteps, with stochastic horizontal diffusion represented by a random-walk term. Land was treated as an impermeable boundary, with particles intersecting land reset to their previous wet position and re-advected with a new random-walk displacement.
Caveats: In the absence of high-resolution hydrodynamic models, a simplified model can provide a structured way to explore likely outcomes, compare scenarios, and guide decisions that would otherwise rely largely on judgement or guesswork. The particle-tracking model below is useful for exploring relative short-term passive dispersal pathways, differences among release sites, retention versus export tendencies, and sensitivity to release timing and horizontal diffusion under the measured logger-depth current conditions, but cannot be used to robustly infer exact particle trajectories, full hydrodynamic circulation, vertical or depth-dependent transport, or wave/wind-driven transport, or small-scale variability in settlement.
Example: Pinnacle Bay, Hook Islands, Whitsundays
Flow fields were generated from TCM-01 tilt logger current observations deployed across the Pinnacle Bay, Hook Islands, Whitsundays in mid-November 2025. Each logger recorded current speed and heading, which were resolved into eastward and northward velocity components and averaged to 60 s intervals. A 12hr window between 11:00 and 23:00 on the 17th November was used to simulate an interpolated flow field.
Code
set.seed(101)
library(tidyverse)
library(janitor)
library(sf)
library(forcats)
library(ggplot2)
tiltmeter_dir <- "/Users/rof011/coralseed/code/tiltmeter"
inst_sf <- tribble(
~lat, ~lon, ~inst,
-20.060709, 148.960138, "2209049",
-20.060966, 148.960434, "2307076",
-20.061180, 148.961006, "2307074",
-20.061194, 148.961482, "2209046",
-20.061212, 148.961916, "1611109",
-20.061279, 148.962260, "2307075",
-20.061236, 148.962533, "2307077",
-20.061206, 148.962770, "2209052",
-20.061131, 148.963162, "2209054",
-20.060709, 148.960138, "22070156",
-20.061194, 148.961482, "1611121",
-20.061206, 148.962770, "1611112"
) |>
mutate(
inst = as.character(inst)
)
tiltmeters <- list.files(
path = tiltmeter_dir,
pattern = "\\.csv$",
full.names = TRUE
) |>
map_dfr(\(file) {
read.csv(file) |>
mutate(
inst = basename(file) |>
str_remove("_.*$")
)
}) |>
janitor::clean_names() |>
mutate(
time = ymd_hms(iso_8601_time, tz = "UTC")
) |>
left_join(inst_sf) |>
drop_na(lon) |>
st_as_sf(
coords = c("lon", "lat"),
crs = 4326,
remove = FALSE
) |>
filter(time >= ymd_hms("2025-11-17T11:00:00.000", tz = "UTC")) |>
filter(time <= ymd_hms("2025-11-17T23:00:00.000", tz = "UTC"))
#
tiltdate <- tiltmeters |> group_by(inst) |> slice(1)
gbr_shp <- dhw::download_gbr_spatial(return="base") |> st_make_valid()
gbr_land <- gbr_shp |> filter(FEAT_NAME %in% c("Island", "Mainland", "Rock")) |> st_transform(4326)
gbr_reef_shp <- gbr_shp |> filter(FEAT_NAME %in% "Reef") |> st_transform(4326)
gbr_buffer <- tiltdate |> st_buffer(250) |> st_bbox()
land_poly <- st_crop(gbr_land, gbr_buffer) |> st_transform(4326)
reef_poly <- st_crop(gbr_reef_shp, gbr_buffer) |> st_transform(4326)
ggplot() + theme_bw() +
geom_sf(data=land_poly, fill="forestgreen") +
geom_sf(data=reef_poly, fill="wheat1", alpha=0.8) +
geom_sf(data=tiltdate, aes(color=inst), show.legend=FALSE) +
geom_sf_text(data=tiltdate, aes(label=inst), size=2.5, position = "jitter", show.legend=FALSE) +
theme(
panel.background = element_rect(fill = "lightblue2", colour = NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(colour = "black", fill = NA, linewidth = 0.5),
plot.margin = margin(0, 0, 0, 0)
) + coord_sf(expand = FALSE) + xlab("Longitude") + ylab("Latitude")
and tiltmeter data in sf format:
Code
iso_8601_time speed_cm_s heading_degrees velocity_n_cm_s
1 2025-11-17T11:00:51.000 0.76 171.89 -0.75
2 2025-11-17T11:01:51.000 0.47 131.96 -0.31
3 2025-11-17T11:02:51.000 0.22 273.45 0.01
4 2025-11-17T11:03:51.000 0.85 143.94 -0.68
5 2025-11-17T11:04:51.000 1.01 142.42 -0.80
6 2025-11-17T11:05:51.000 1.74 132.94 -1.19
velocity_e_cm_s inst time lat lon
1 0.11 1611109 2025-11-17 11:00:51 -20.06121 148.9619
2 0.35 1611109 2025-11-17 11:01:51 -20.06121 148.9619
3 -0.22 1611109 2025-11-17 11:02:51 -20.06121 148.9619
4 0.50 1611109 2025-11-17 11:03:51 -20.06121 148.9619
5 0.61 1611109 2025-11-17 11:04:51 -20.06121 148.9619
6 1.28 1611109 2025-11-17 11:05:51 -20.06121 148.9619
geometry
1 148.96192, -20.06121
2 148.96192, -20.06121
3 148.96192, -20.06121
4 148.96192, -20.06121
5 148.96192, -20.06121
6 148.96192, -20.06121
flow_speed across the sample window (tilt-meters ordered west-east across the bay)
Interpolate flow field
The simulate_flowfields() function takes point-based velocity observations from tiltmeters and spatially interpolates across the study domain using inverse distance weighting (gstat::idw()) to create a time-varying two-dimensional flow field (user defined, 10m spatial resolution) at each timestep (user defined, 5 minute resolution). Interpolation was restricted to the ocean domain, with land polygons used to define barriers and exclude non-ocean areas. The resulting gridded flow fields describe the interpolated 2D current structure across the logger array. As the flow fields are derived from spatial interpolation of fixed-point measurements rather than a hydrodynamic model, they represent an observation-driven approximation of near-bed or logger-depth currents during the deployment period.
Code
flow_sim <- coralseed::simulate_flowfields(
tiltmeters = tiltmeters,
gbr_buffer = gbr_buffer,
crs_m = 32755,
grid_res_m = 10,
buffer_m = 250,
time_unit = "5 minutes",
idp = 2
)
Example tilt-meter derived flow-field at 11:30:
Code
library(terra)
grid_res <- 10
r_template <- terra::rast(
terra::ext(flow_sim$grid_sf),
resolution = grid_res,
crs = terra::crs(flow_sim$grid_sf)
)
flow_raster <- flow_sim$flow_field |>
filter(time == ymd_hms("2025-11-17T11:30:00", tz = "UTC")) |>
mutate(
xend = x + u_ms * 10,
yend = y + v_ms * 10
) |>
vect() |>
select(speed_ms) |>
rasterize(y=r_template, field = "speed_ms") |>
project(crs(land_poly)) #|>
#crop(land_poly, mask=FALSE)
flow_plot <- flow_sim$flow_field |>
filter(time == ymd_hms("2025-11-17T11:30:00", tz = "UTC")) |>
mutate(
xend = x + u_ms * 1200,
yend = y + v_ms * 1200
)
ggplot() +
tidyterra::geom_spatraster(
data = flow_raster,
alpha=0.8
) +
geom_segment(
data = flow_plot |> st_drop_geometry(),
aes(x = x, y = y, xend = xend, yend = yend),
arrow = arrow(length = unit(0.03, "inches")),
linewidth = 0.15,
alpha = 0.8,
colour = "black"
) +
geom_sf(data = st_transform(land_poly, 32755), fill = "beige") +
coord_sf(crs = st_crs(32755)) +
theme_bw() +
xlab("Longitude") + ylab("Latitude") +
scale_fill_distiller(palette="RdBu", na.value = "transparent", name = "speed (m/s)")
Seed particles
simulate_particles() simulates passive particle dispersal using a two-dimensional Lagrangian particle-tracking approach driven by interpolated current-meter observations. Eastward and northward velocity components from the tiltmeters are averaged to 60 s intervals and spatially interpolated across the bay to generate a time-varying flow field. Particles are released from user-defined fixed source locations (here n=3) at a user-defined rate (here n=1,000 particles) at a fixed rate (here every 60 s for 1 h, n=180,000 particles in total). Particle positions were updated at 60 s timesteps using local interpolated velocities, with horizontal diffusion represented by a random-walk term, where *K** is the horizontal diffusion coefficient integrated into a random displacement per timestep (rnorm(n(), mean = 0, sd = sqrt(2 * K * dt))). A random-walk displacement was applied independently in the x and y directions at each timestep, with default standard deviation calculated from the horizontal diffusivity coefficient as sqrt(2 * K * dt) (i.e. particle movemement per timestep = current displacement + random diffusion). For the simulations reported here, K = 0.05 m² s⁻¹ and dt = 60 s, giving a per-timestep random-walk standard deviation of sqrt(2 * 0.05 * 60) = 2.45 m. The model also allows this displacement scale to be specified directly using random_walk_sd.
To avoid “beached” particles, when a particle intersects with landits position is reset to the previous “wet” location allowing it to be re-advected on the next timestep while retaining the stochastic diffusion term. This iterative re-release/random-walk procedure prevents particles from becoming trapped on land boundaries while maintaining land as an impermeable barrier to net transport. Alternatively, this can be switched of and beached particles are removed from the model. Importantly, the model represents a simplified, observation-driven 2D Lagrangian advection–diffusion simulation rather than a full hydrodynamic model.
*K can be set from K = 0, where particles follow only the interpolated currents with no diffusion, through weak (K = 0.01), modest small-bay (K = 0.05), and stronger (K = 0.1) spreading, up to very strong coastal mixing (K = 1).
Note - simulate_particles() uses furrr and future to simulate independent particle groups in parallel as needed (180k individual simulated particles runs in ~2.4mins)
Code
# set release sites (here n=3)
release_sites <- tribble(
~lat, ~lon,
-20.0612, 148.9605,
-20.0614, 148.9615,
-20.06125, 148.9625
) |>
st_as_sf(
coords = c("lon", "lat"),
crs = 4326,
remove = FALSE
)
# simulate particles
particle_sim <- coralseed::simulate_particles(
flow_field = flow_sim$flow_field,
land_poly = land_poly,
release_sites = release_sites,
crs_m = 32755,
dt = 5 * 60,
release_duration_seconds = 60 * 60,
n_particles_per_site_per_release = 100,
K = 0.05,
max_land_retry = 200,
seed = 101,
parallel = TRUE,
workers = 6
)
Below are n=12 individual subsampled particle tracks from the particle dispersal model
Code
track_sample <- particle_sim$particle_tracks_sf |>
filter(particle_id %in% sample(unique(particle_id), 12))
ggplot() +
theme_bw() +
theme(
panel.background = element_rect(fill = "lightblue2", colour = NA),
#plot.background = element_rect(fill = "lightblue2", colour = NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(colour = "black", fill = NA, linewidth = 0.5),
plot.margin = margin(0, 0, 0, 0)
) +
geom_sf(
data = land_poly,
fill = "forestgreen",
colour = NA
) +
geom_sf(data=reef_poly, fill="wheat1", alpha=0.8) +
geom_path(
data = track_sample |> st_drop_geometry(),
aes(
x = x,
y = y,
group = particle_uid
),
color="black",
alpha = 1,
linewidth = 1
) +
geom_sf(
data = track_sample |> slice_max(n=1, order_by=time),
aes(group = particle_uid),
color="black", shape=21,
size = 2) +
geom_sf(
data = track_sample |> slice_min(n=1, order_by=time),
aes(group = particle_uid),
color="black", shape=21, color="black",
size = 2) +
geom_sf(
data = release_sites,
color="black", stroke=0.2,
size = 1, alpha=0.2
) +
# scale_color_distiller(palette="Spectral") +
coord_sf(crs = st_crs(track_sample), expand = FALSE) +
labs(
colour = "Release site",
x = NULL,
y = NULL
) +
facet_wrap(~particle_uid)
Combined particle tracks for a subset of 100 particle tracks for the full 12 hour period highlighting i) local retention within the bay across the 12hr period, ii) advection of currents and loss of a limited number of particles from the western point, and iii) loss of particles through the eastern passage returning through the rocks through a change in tidal currents ~8hrs post-release
Code
land_poly_2 <- st_crop(gbr_land |> st_transform(st_crs(particle_sim$particle_tracks_sf)),
particle_sim$particle_tracks_sf |> st_buffer(100) |> st_bbox()) |> st_transform(4326)
reef_poly_2 <- st_crop(gbr_reef_shp |> st_transform(st_crs(particle_sim$particle_tracks_sf)),
particle_sim$particle_tracks_sf |> st_buffer(100) |> st_bbox()) |> st_transform(4326)
track_sample2 <- particle_sim$particle_tracks_sf |>
filter(particle_id %in% sample(unique(particle_id), 100))
ggplot() +
theme_bw() +
theme(
panel.background = element_rect(fill = "lightblue2", colour = NA),
#plot.background = element_rect(fill = "lightblue2", colour = NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(colour = "black", fill = NA, linewidth = 0.5),
plot.margin = margin(0, 0, 0, 0)
) +
geom_sf(
data = land_poly_2,
fill = "forestgreen",
colour = NA
) +
geom_sf(data=reef_poly, fill="wheat1", alpha=0.8) +
geom_path(
data = track_sample2 |> st_drop_geometry(),
aes(
x = x,
y = y,
group=particle_uid,
color = age_seconds/60
),
alpha = 1,
linewidth = 0.1
) +
geom_sf(
data = release_sites,
color="black", #stroke=0.2,
size = 2, alpha=0.5
) +
# scale_color_distiller(palette="RdBu") +
scale_color_viridis_c(option = "magma", direction=-1) +
coord_sf(crs = st_crs(particle_sim$particle_tracks_sf), expand = FALSE) +
labs(
colour = "Time\n (mins)",
x = NULL,
y = NULL
)
The model can be run with or without “beaching”: rerunning the above simulation results in >99% of particles becoming beached and removed from the simulation due to high numbers of particles hitting land in the high-circulation recirculating environment of Pinnacle Bay.
Code
# simulate particles
particle_sim_beached <- coralseed::simulate_particles(
flow_field = flow_sim$flow_field,
land_poly = land_poly,
release_sites = release_sites,
crs_m = 32755,
dt = 5 * 60,
release_duration_seconds = 60 * 60,
n_particles_per_site_per_release = 100,
K = 0.05,
max_land_retry = 0,
seed = 101,
parallel = TRUE,
workers = 6,
beached = TRUE
)
particle_sim_beached$particle_track_summary |> slice(1-5)
# A tibble: 1 × 16
n_particles n_particle_uids n_track_rows n_beached n_beached_particles
<int> <int> <int> <int> <int>
1 3600 3600 122398 3599 3599
# ℹ 11 more variables: prop_beached_particles <dbl>, mean_land_retry_n <dbl>,
# max_land_retry_n <int>, prop_retried_land <dbl>, dt <dbl>, K <dbl>,
# random_walk_sd <dbl>, random_walk_sd_source <chr>, parallel <lgl>,
# workers <dbl>, beached <lgl>
