In pbs-assess/PACea: An R package of Pacific ecosystem information to help facilitate an ecosystem approach to fisheries management

rmarkdown::render("hotssea.Rmd")
# to build, or click the knit button in RStudio
# Can take around 10 minutes because of the big panel figure

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#",
  cache = TRUE,
  cache_path = "hotssea-cache/",
  fig.path = "hotssea-figs-cache/",
  fig.width = 5.7,
  fig.height = 5
)

Introduction

The Hindcast of the Salish Sea (HOTSSea) is a physical ocean model that recreates conditions throughout the Salish Sea from 1980 to 2018, filling in the gaps in patchy measurements Oldford et al., (in review). The model predicts physical ocean properties with sufficient accuracy to be useful for a variety of applications. It corroborates observed ocean temperature trends and can examine areas with few observations. Results indicate that some seasons and areas are warming faster than others. HOTSSea was developed using the Nucleus for European Modelling of the Ocean (NEMO) engine, and temperature and salinty outputs are included in pacea (from version 1 of HOTSSea).

The HOTSSea outputs in pacea are adapted from results provided by Greig Oldford (Fisheries and Oceans Canada). For model details see Oldford et al., (in review), which should be read and cited if you use the results, and Greig and Andrew Edwards should be consulted for any clarifications and uses of the results.

The model domain includes several distinct geographic areas within the Salish Sea: the Juan de Fuca Strait, Strait of Georgia, Gulf Islands, and Puget Sound. The model grid is approximately 1.5 km in horizontal resolution and has a width of ~200 km and length of ~450 km (132 cells x 299 cells) and was rotated 29° counter clockwise to true north to align with the axis of the Strait of Georgia. For pacea, results were translated onto a north-south aligned grid of 2 km resolution (all such calculations are in https://github.com/pbs-assess/data-raw/hotssea/hotssea-data-interpolation.R), resulting in 6165 spatial cells. The vertical grid for the HOTSSea model is divided into 40 vertical levels that increase in resolution at the surface. Monthly averages of variables over several specified depth ranges are included in pacea, for every month from January 1980 to December 2018.

We have wrangled the HOTSSea model outputs into the same format as those from the British Columbia Continential Margin (BCCM) model, such that the functions regarding plotting, climatology, and anomalies, described in the BCCM vignette will work for the HOTSea model output. Therefore, in this vignette we describe the available HOTSSea model outputs and present some results, and refer users to see the BCCM vignette for further analytical and visual possibilities (to avoid duplicating details here).

Furthermore, the grid for the HOTSSea outputs matches that used for the full BCCM outputs; see the BCCM_full vignette.

library(pacea)
library(dplyr)
library(sf)
library(ggplot2)

Available HOTSSea results

As for the BCCM results, the HOTSSea results were wrangled to simplify data format and structure. Because the results are spatially and temporally resolved, the compressed file size of each variable is fairly large (around 8 Mb, and there are 40 files). To keep the pacea package a reasonable downloadable size, the HOTSSea results are stored on Zenodo and can be downloaded locally using functions within pacea (obviously requiring internet access). The downloading is only needed once as the files get cached locally on your computer.

The list of variables of the HOTSSea results available can be viewed from the list hotssea_data:

hotssea_data

Note that the naming convention of variables is the same as for the BCCM results, such that avg0to30m means the average over all depths from 0 to 30 m, for example. The suffixes _min, _mean, _max, and _std are additionally used here.

As an explicit example, hotssea_avg0to30m_temperature_min is calculated as follows. For each cell in the HOTSSea model three-dimensonal domain:

• the HOTSSea model calculates temperature on a time step of two minutes
• the model saves the results as hourly means of the two-minute values
• the daily mean of those hourly means is calculated
• the minimum daily mean over the given month is calculated.

Then:

• for calculations that are averages over a depth range, such as hotssea_avg0to30m_temperature_min, for each (two-dimensional) spatial location the variable is average of the values over all cells within that depth range.
• for surface or bottom variables, such as hotssea_bottom_temperature_min, the variable is the value in the surface or bottom cell for each (two-dimensional) spatial location.

Similar calculations are done for the _mean (the mean daily mean), _max (the maximum of the daily means), and _std (the standard deviation of the daily means), replacing the word 'minimum' in the fourth bullet point. And the same calculations are done for salinity.

So, hotssea_avg0to30m_temperature_min is the depth-integrated mean (over the top 30 m) of each modelled depths' minimum (over the month) daily mean temperature. And hotssea_bottom_salinity_sd is the the bottom cell's standard deviation (over the month) of daily mean salinity.

The way the averaging was done means that the within-day variability of temperature and salinity are not retained (and averaging over depth ranges also sacrifices some variability). The alternative would be to base the monthly statistics on the hourly model outputs (skipping the third step above). Both options have strengths and drawbacks. Variability on hourly time scales may be important for some ecological groups, possibly zooplankton, but on the other hand may result in monthly fields that are too noisy for research focused on other ecological groups, like fish. So we opted for the given approach. With some work the within-day variability quantities could be calculated (contact us).

Model depths vary by grid cell, usually only slightly at the surface. So variables are therefore processed accounting for varying depth levels. See the help file for any HOTSSea variable (e.g. ?hotssea_surface_salinity_min) for explicit details of depth resolution, including how surface values were calculated.

Downloading all HOTSSea results

It is simplest to do a one-time download of all HOTSSea results. This will take several minutes (and uses parallel processing if your computer has the capability). Simply run:

hotssea_all_variables()

This downloads from our Zenodo archive all variables into a cache folder. This is identified using:

pacea_cache()
"C:\\Users\\<your-windows-username>\\AppData\\Local/pacea/Cache"
# On linux and Macs this will look different.

and the HOTSSea results will be placed in the hotssea subfolder. See ?hotssea_all_variables if you have problems.

Using the HOTSSea results

For example, the hotssea_avg0to30m_temperature_max() function gives the HOTSSea model's mean (over the top 30 m) of each modelled depths' maximum (over the month) daily mean temperature:

hot <- hotssea_avg0to30m_temperature_max()

Note the (). The hotssea_avg0to30m_temperature_max() function assigns the output to our designated variable simply named hot.

If you have not already previously downloaded the results it will first download them into your cache folder (i.e. it will get just the specified model results, rather than all of them as in the hotssea_all_variables() described above).

View the data help for more information on HOTSSea results (e.g. ?hotssea_avg0to30m_temperature_max); these all point to a common help file for all hotssea objects.

So, what does the object look like?

hot
head(hot[, 1:5])   # note that the geometry column is included when selecting columns from
                   # an `sf` object.

The object is in wide format, with each column representing a unique year-month combination, and is a 'simple features' (sf) R object, and so contains a geometry column. The data also have various attributes, such as the units for the data values, and some extra ones for HOTSSea used to automate some plotting.

names(attributes(hot))
attributes(hot)$units
class(hot)

HOTSSea visualization

HOTSSea results can be plotted using the plot() function for a quick view of the data. By using plot(), since the results object has our class pacea_st (same as for BCCM), plot() automatically uses our customised plot.pacea_st(). The default settings are to plot April 2018:

plot(hot)

The code works the same as in the BCCM vignette, for example showing two months for each of three years: and

plot(hot,
     months = c("January", "September"),
     years = c(1995, 2010, 2018),
     eez = FALSE)         # Do not show the dotted line for the EEZ

To inspect the full time series we can show every available month in a huge panel plot (which takes a while to plot). We are looking at hotssea_avg0to30m_temperature_max (from the hot <- command above). So the following figure shows, for each month from 1980 to 2018, the depth-integrated mean(over the top 30 m) of each modelled depths' maximum (over the month) daily mean temperature.

plot(hot, months = 1:12,
     years = 1980:2018,
     eez = FALSE)

Calculating climatogy and anomalies of HOTSSea results

The pacea package is about more than producing pretty pictures, since all the data and model output is available for your own analyses.

The package has some built in functions to calculate climatologies and anomalies of the spatiotemporal results, and these functions can be used on the HOTSSea results. See help files for details of functions. Here are some brief examples taken from the BCCM vignette.

So let's ask: how do the anomalies of minimum bottom temperature for April and September in 2010 and 2018, look compared to a climatology (for each month) from 1980 to 2000. The minimum bottom temperature is defined specifically as the bottom cell's minimum (over the month) daily mean salinity.

First, let's load the model results and calculate the climatology, which returns data in long format.

hot_bottom <- hotssea_bottom_temperature_min()
clim <- calc_clim(hot_bottom,
                  clim_years = 1980:2000)
clim   # Each row is a month-polygon combination

The anomalies relative to a climatology can be calculated in one function:

anom <- calc_anom(hot_bottom,
                  clim_years = 1980:2000,
                  time_period_return = c("Apr", "Sep"),
                  years_return = c(2010, 2018))
anom
plot(anom,
     months.plot = c("Apr", "Sep"),
     years.plot = c(2010, 2018),
     bc = FALSE,
     eez = FALSE)

We have also incorporated the ability to add climatological data to the plot, which creates contour lines showing the anomaly areas that are above (or below for other variables) the 90th and 99th percentile of data.

See the rest of the BCCM vignette for details and other features, like estimating statistics across spatial areas.

Updating HOTSSea data

When you download HOTSSea data, the most recent version of the data will be downloaded to your cache folder. Occasionally, there might be updates to the HOTSSea data (e.g. adding data from more recent years and maybe future projections). If you want to check for an update for a data file you have already downloaded to the pacea_cache directory: Keep an eye on the NEWS file on GitHub for updates (for such major updates will be also likely mention them on the main README). For now (November 2024) there are no immediate plans for updates.

Reference

Oldford, G.L., Jarníková, T., Christensen, V., and Dunphy, M. (in review). HOTSSea v1: a NEMO-based physical Hindcast of the Salish Sea (1980–2018) supporting ecosystem model development. Preprint. https://doi.org/10.5194/gmd-2024-58 .

pbs-assess/PACea documentation built on April 17, 2025, 11:36 p.m.