4.9 Wave benthic orbital velocities pt 2
4.9 Wave benthic orbital velocities pt 2
On this page
This page is part 2 of a two-part wave series. It reads the significant wave height (Hsig) and wave period (Per) series that part 1 (4.8) extracted, then converts them into benthic orbital velocities (BOVs) for each of the 50 USVI reef sites. The BOV calculation solves the linear wave dispersion relation for wavelength at each site depth, applies the Smith et al. 2016 orbital-velocity equation (Smith et al. 2016), and runs across cores in parallel. It produces per-site BOV summary statistics (mean, percentiles, decile mean) and the raw per-timestep BOV series, both offered for download below and used as external drivers in downstream resilience work.
Data sources
This page reads the per-coordinate Hsig/Per series and the site-to-grid-coordinate lookup written by part 1, which derive from the waves dataset (SWAN aggregation model), together with a coastline shapefile used for the land check. The wave record spans January 1, 2015 through December 31, 2020. The per-site BOV summary and raw BOV series produced here are available in the Downloads section below.
Import data
We load the master catalog of 50 USVI reef sites into sitedat. This catalog comes from the waves subfolder because it carries the grid coordinate we determined for each site in an earlier version of this analysis.
We read the Hsig and Per series that part 1 wrote, drop duplicate time-by-coordinate rows, parse the timestamps, and sort by time. We also load the site-to-coordinate lookup.
Show code
data <- read.csv("../../outputs/s4pt8_SWANHsigPer_AllSiteCoords_20150101_20201231.csv")
#remove duplicate pairs of time and site_coord
data <- data[!duplicated(data[,c("time","site_coord")]),]
data$time <- as.POSIXct(data$time, format = "%m/%d/%Y %H:%M:%S")
data <- data |> arrange(time)
sitecoords <- read.csv("../../outputs/s4pt8_SWANHsigPer_sitecoords.csv")The calcbov function
What happens here: calcbov() computes benthic orbital velocities (bov) from wave Hsig and Per, following the equations in Smith et al. 2016 (Smith et al. 2016).
Setup. The bovraw data frame stores the benthic orbital velocities (BOV) across sites and times. Gravitational acceleration g is 9.81 m/s².
Functionality.
Read slices of significant wave height (
Hsig) and wave period (Per) data.Solve the following linear dispersion relation for L.
\[ \omega^2 = gk \tanh(kh)\]
Here \(\omega\) is wave angular frequency, \(g\) is acceleration due to gravity, \(k\) is wavenumber, and \(h\) is depth.
We substitute \(\omega\) and \(k\) with their equivalent expressions in terms of wave period \(T_p\) (\(\omega = 2\pi/T_p\)) and wavelength \(L\) (\(k = 2\pi/L\)):
\[(\frac{2\pi}{T})^2 = g(\frac{2\pi}{L}) \tanh(\frac{2\pi}{L}h)\]
Solving this equation means finding its root: the value of \(L\) for which the two sides are equal. We define a function \(f(L)\) that computes the left-hand side minus the right-hand side:
\[f(L) = (\frac{2\pi}{T})^2 - g(\frac{2\pi}{L}) \tanh(\frac{2\pi}{L}h) \]
The root of \(f(L)\) (the value of \(L\) that makes \(f(L) = 0\)) satisfies the original equation.
-
Calculate BOV ( \(u_b\) ) using:
\[u = \frac{H g T_p}{2 L \cosh\left(\frac{2 \pi d}{L}\right)}\]
Output.
bov
Parallelization. We set up a parallel backend with the parallel package to speed up the computation. calcbov then runs in parallel across cores, which cuts the total computation time substantially (last measured at 8 s for parallel execution against 34 s for sequential execution).
Show code
# Initialize bovraw
bovraw <- data.frame(matrix(ncol = nrow(sitedat), nrow = length(unique(data$time))))
colnames(bovraw) <- sitedat$site
rownames(bovraw) <- unique(data$time)
g <- 9.81
# Function for solving for L with retry mechanism
solve_L <- function(T, g, h) {
f <- function(L) {
((2 * pi) / T) ^ 2 - g * ((2 * pi) / L) * tanh(((2 * pi) / L) * h)
}
intervals <- list(
c(0.01, 100),
c(0.1, 100),
c(1, 100),
c(0.01, 200),
c(0.1, 200),
c(1, 200)
)
for (interval in intervals) {
solution <- tryCatch({
uniroot(f, interval = interval, tol = 1e-6)
}, error = function(e) {
message("Error in uniroot with interval ", paste(interval, collapse = ", "), ": ", e)
return(NULL)
})
if (!is.null(solution)) {
return(solution$root)
}
}
message("Failed to find root for T = ", T, ", h = ", h)
return(NA)
}
# Function to calculate bov
calcbov <- function(i, sitecoords, sitedat, data, g) {
sitind <- which(sitedat$site == colnames(bovraw)[i])
s <- sitecoords$site[sitind]
sc <- sitecoords$site_coord[sitind]
d <- sitedat$depth[which(sitedat$site == s)]
dt <- data[which(data$site_coord == sc), ]
hs <- dt$hsig
p <- dt$per
t <- dt$time
L_wavelength <- sapply(p, function(x) solve_L(x, g, d))
bov <- (hs * g * p) / (2 * L_wavelength * cosh(2 * pi * d / L_wavelength))
bov <- data.frame(bov)
colnames(bov) <- s
return(bov)
}
# Measure time for parallel execution
parallel_time <- system.time({
numCores <- detectCores() - 1
cl <- makeCluster(numCores)
clusterExport(cl, c("sitecoords", "sitedat", "data", "g", "solve_L", "calcbov", "bovraw"))
results <- parLapply(cl, 1:ncol(bovraw), function(i) {
calcbov(i, sitecoords, sitedat, data, g)
})
combined_results <- do.call(cbind, lapply(results, function(df) df[[1]]))
site_names <- sapply(results, function(df) colnames(df))
colnames(combined_results) <- site_names
rownames(combined_results) <- rownames(bovraw)
bovraw <- combined_results
stopCluster(cl)
})To compare against the sequential run time, uncomment the following code block.
Show code
# Measure time for sequential execution
# sequential_time <- system.time({
# sequential_bovraw <- data.frame(matrix(ncol = nrow(sitedat), nrow = length(unique(data$time))))
# colnames(sequential_bovraw) <- sitedat$site
# rownames(sequential_bovraw) <- unique(data$time)
#
# results <- lapply(1:ncol(sequential_bovraw), function(i) {
# calcbov(i, sitecoords, sitedat, data, g)
# })
# combined_results <- do.call(cbind, lapply(results, function(df) df[[1]]))
# site_names <- sapply(results, function(df) colnames(df))
# colnames(combined_results) <- site_names
# rownames(combined_results) <- rownames(sequential_bovraw)
# sequential_bovraw <- combined_results
# })
#
# print(sequential_time)Summarize BOV per site
For each site we reduce the per-timestep BOV series into summary statistics: the mean, the 5th, 50th, and 95th percentiles, the absolute minimum and maximum, and the mean of the top decile. We attach each site’s depth and coordinates so the summary can be joined and mapped downstream.
Show code
summary_stats <- lapply(colnames(bovraw), function(site) {
bovs <- bovraw[,which(colnames(bovraw) == site)]
depth <- sitedat$depth[which(sitedat$site == site)]
lat <- sitecoords$lat[which(sitecoords$site == site)]
lon <- sitecoords$lon[which(sitecoords$site == site)]
gridLat <- sitecoords$gridLat[which(sitecoords$site == site)]
gridLon <- sitecoords$gridLon[which(sitecoords$site == site)]
newRow <- data.frame(
site = site,
depth = depth,
lat = lat,
lon = lon,
gridLat = gridLat,
gridLon = gridLon,
bov_average = mean(bovs, na.rm = TRUE),
bov_5percentile = quantile(bovs, 0.05, na.rm = TRUE),
bov_50percentile = median(bovs, na.rm = TRUE),
bov_95percentile = quantile(bovs, 0.95, na.rm = TRUE),
bov_max = max(bovs, na.rm = TRUE),
bov_min = min(bovs, na.rm = TRUE),
top_decile_mean = mean(bovs[order(-bovs)][1:floor(0.10 * length(bovs))], na.rm = TRUE)
)
return(newRow)
})
# Combine the list into a single data frame
summary_stats_df <- do.call(rbind, summary_stats)Check sites with NA BOVs
Some sites can return NA BOVs when their grid point falls on land. We check that here by mapping the grid points and comparing them to the coastline.
We load the coastline shapefile from the USGS basemaps and crop it to a bounding box around the USVI.
Show code
#bounding box coordinates
xmin <- -65.1 # min(sitedat$lon)
ymin <- 17.6 # min(sitedat$lat)
xmax <- -64.4 # max(sitedat$lon)
ymax <- 18.4 # max(sitedat$lat)
# import the shapefile as an sf object
my_shapefile_sf <- st_read("../shapefiles/pvishrpl/pvishrpl.shp", quiet = TRUE)
# Crop the sf object
cropped_sf <- st_crop(my_shapefile_sf, xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax)
# the grid itself We pull the gridded coordinates for the sites with NA BOV values so we can check their distance to land.
Show code
bov_site <- summary_stats_df
bov_site_NA <- bov_site[which(is.na(bov_site$bov_average)), ]
ptsok <- data.frame(x1 = bov_site$gridLon, y1 = bov_site$gridLat)
sptsok <- SpatialPoints(list(x = ptsok$x1, y = ptsok$y1))
pts <- data.frame(x1 = bov_site_NA$gridLon, y1 = bov_site_NA$gridLat)
spts <- SpatialPoints(list(x = pts$x1, y = pts$y1))The map below plots every grid point over the USVI coastline so we can see which ones sit on land. Green points return a valid BOV; red points return NA, and a red point over land confirms the grid cell fell on the coast rather than in the water. This is a spatial diagnostic, so it is built as an interactive ggplotly map that supports panning and zooming. The chunk is prepared but not evaluated in the published render (eval: false); enable it locally to inspect the land overlaps.
Show code
# Grid points for the 50 USVI reef sites over the coastline. Green points return
# a valid benthic orbital velocity; red points return NA, which flags a grid cell
# that fell on land rather than in the water. Pan and zoom to inspect sites.
# plot
myplot <- ggplot() +
geom_sf(data = cropped_sf, fill = "grey", color = "grey", size = 0.25) +
geom_point(data=ptsok, aes(x=x1, y=y1), color = "green", size = 0.2) +
geom_point(data=pts, aes(x=x1, y=y1), color = "red", size = 0.2) +
theme_void() # This removes axes, similar to par(mar = c(0, 0, 0, 0))
# Convert to an interactive plotly plot
ggplotly(myplot)Downloads
The links below download the per-site BOV summary, the raw per-timestep BOV series, and the metadata for each, all generated from the objects built above.
- Derived data:
s4pt9_bov_summary_statistics_50sites.csv - Metadata:
s4pt9_bov_summary_statistics_50sites.txt - Derived data:
s4pt9_bov_raw_50sites_2015_2020.csv - Metadata:
s4pt9_bov_raw_50sites_2015_2020.txt
version 1.0.0 • in-review • data ≤ 2023-12-31