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.
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
# data_name
# 1 hotssea_surface_salinity_min
# 2 hotssea_surface_salinity_mean
# 3 hotssea_surface_salinity_max
# 4 hotssea_surface_salinity_std
# 5 hotssea_surface_temperature_min
# 6 hotssea_surface_temperature_mean
# 7 hotssea_surface_temperature_max
# 8 hotssea_surface_temperature_std
# 9 hotssea_avg0to30m_salinity_min
# 10 hotssea_avg0to30m_salinity_mean
# 11 hotssea_avg0to30m_salinity_max
# 12 hotssea_avg0to30m_salinity_std
# 13 hotssea_avg0to30m_temperature_min
# 14 hotssea_avg0to30m_temperature_mean
# 15 hotssea_avg0to30m_temperature_max
# 16 hotssea_avg0to30m_temperature_std
# 17 hotssea_avg30to150m_salinity_min
# 18 hotssea_avg30to150m_salinity_mean
# 19 hotssea_avg30to150m_salinity_max
# 20 hotssea_avg30to150m_salinity_std
# 21 hotssea_avg30to150m_temperature_min
# 22 hotssea_avg30to150m_temperature_mean
# 23 hotssea_avg30to150m_temperature_max
# 24 hotssea_avg30to150m_temperature_std
# 25 hotssea_avg150toBot_salinity_min
# 26 hotssea_avg150toBot_salinity_mean
# 27 hotssea_avg150toBot_salinity_max
# 28 hotssea_avg150toBot_salinity_std
# 29 hotssea_avg150toBot_temperature_min
# 30 hotssea_avg150toBot_temperature_mean
# 31 hotssea_avg150toBot_temperature_max
# 32 hotssea_avg150toBot_temperature_std
# 33 hotssea_bottom_salinity_min
# 34 hotssea_bottom_salinity_mean
# 35 hotssea_bottom_salinity_max
# 36 hotssea_bottom_salinity_std
# 37 hotssea_bottom_temperature_min
# 38 hotssea_bottom_temperature_mean
# 39 hotssea_bottom_temperature_max
# 40 hotssea_bottom_temperature_stdNote 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:
Then:
hotssea_avg0to30m_temperature_min, for each
(two-dimensional) spatial location the variable is average of the values
over all cells within that depth range.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.
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:
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.
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:
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
# Simple feature collection with 6165 features and 468 fields
# Geometry type: POLYGON
# Dimension: XY
# Bounding box: xmin: 988882.9 ymin: 232740 xmax: 1284883 ymax: 656740
# Projected CRS: NAD83 / BC Albers
# # A tibble: 6,165 × 469
# `1980_1` `1980_2` `1980_3` `1980_4` `1980_5` `1980_6` `1980_7` `1980_8`
# * <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 10.2 9.82 9.97 10.2 11.4 11.9 12.6 12.1
# 2 10.2 9.85 10.0 10.3 11.5 12.1 12.9 12.3
# 3 10.2 9.82 9.96 10.2 11.3 11.9 12.6 12.1
# 4 10.2 9.84 9.98 10.2 11.3 11.9 12.6 12.1
# 5 10.3 9.89 10.1 10.3 11.5 12.1 12.9 12.3
# 6 10.2 9.83 9.95 10.2 11.2 11.9 12.6 12.1
# 7 10.2 9.84 9.96 10.2 11.2 11.9 12.6 12.0
# 8 10.3 9.90 10.0 10.2 11.3 11.9 12.6 12.0
# 9 10.2 9.84 9.93 10.2 11.2 11.9 12.6 12.0
# 10 10.3 9.89 9.96 10.2 11.2 11.9 12.6 12.0
# # ℹ 6,155 more rows
# # ℹ 461 more variables: `1980_9` <dbl>, `1980_10` <dbl>, `1980_11` <dbl>,
# # `1980_12` <dbl>, `1981_1` <dbl>, `1981_2` <dbl>, `1981_3` <dbl>,
# # `1981_4` <dbl>, `1981_5` <dbl>, `1981_6` <dbl>, `1981_7` <dbl>,
# # `1981_8` <dbl>, `1981_9` <dbl>, `1981_10` <dbl>, `1981_11` <dbl>,
# # `1981_12` <dbl>, `1982_1` <dbl>, `1982_2` <dbl>, `1982_3` <dbl>,
# # `1982_4` <dbl>, `1982_5` <dbl>, `1982_6` <dbl>, `1982_7` <dbl>, …
head(hot[, 1:5]) # note that the geometry column is included when selecting columns from
# Simple feature collection with 6 features and 5 fields
# Geometry type: POLYGON
# Dimension: XY
# Bounding box: xmin: 1076883 ymin: 650740 xmax: 1084883 ymax: 656740
# Projected CRS: NAD83 / BC Albers
# # A tibble: 6 × 6
# `1980_1` `1980_2` `1980_3` `1980_4` `1980_5` geometry
# <dbl> <dbl> <dbl> <dbl> <dbl> <POLYGON [m]>
# 1 10.2 9.82 9.97 10.2 11.4 ((1080883 656740, 1082883 656740…
# 2 10.2 9.85 10.0 10.3 11.5 ((1082883 656740, 1084883 656740…
# 3 10.2 9.82 9.96 10.2 11.3 ((1078883 654740, 1080883 654740…
# 4 10.2 9.84 9.98 10.2 11.3 ((1080883 654740, 1082883 654740…
# 5 10.3 9.89 10.1 10.3 11.5 ((1082883 654740, 1084883 654740…
# 6 10.2 9.83 9.95 10.2 11.2 ((1076883 652740, 1078883 652740…
# 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.
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:
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 EEZTo 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.
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
# Simple feature collection with 73980 features and 4 fields
# Geometry type: POLYGON
# Dimension: XY
# Bounding box: xmin: 988882.9 ymin: 232740 xmax: 1284883 ymax: 656740
# Projected CRS: NAD83 / BC Albers
# # A tibble: 73,980 × 5
# month clim_value clim_sd clim_n geometry
# <dbl> <dbl> <dbl> <int> <POLYGON [m]>
# 1 1 8.95 0.191 21 ((1080883 656740, 1082883 656740, 1082883 65…
# 2 2 8.94 0.187 21 ((1080883 656740, 1082883 656740, 1082883 65…
# 3 3 8.94 0.186 21 ((1080883 656740, 1082883 656740, 1082883 65…
# 4 4 8.94 0.188 21 ((1080883 656740, 1082883 656740, 1082883 65…
# 5 5 8.94 0.190 21 ((1080883 656740, 1082883 656740, 1082883 65…
# 6 6 8.94 0.193 21 ((1080883 656740, 1082883 656740, 1082883 65…
# 7 7 8.93 0.195 21 ((1080883 656740, 1082883 656740, 1082883 65…
# 8 8 8.92 0.193 21 ((1080883 656740, 1082883 656740, 1082883 65…
# 9 9 8.93 0.191 21 ((1080883 656740, 1082883 656740, 1082883 65…
# 10 10 8.93 0.188 21 ((1080883 656740, 1082883 656740, 1082883 65…
# # ℹ 73,970 more rowsThe 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
# Simple feature collection with 6165 features and 4 fields
# Geometry type: POLYGON
# Dimension: XY
# Bounding box: xmin: 988882.9 ymin: 232740 xmax: 1284883 ymax: 656740
# Projected CRS: NAD83 / BC Albers
# # A tibble: 6,165 × 5
# `2010_4` `2010_9` `2018_4` `2018_9` geometry
# <dbl> <dbl> <dbl> <dbl> <POLYGON [m]>
# 1 0.00324 -0.0210 0.280 0.300 ((1080883 656740, 1082883 656740, 108288…
# 2 0.244 0.138 0.336 0.349 ((1082883 656740, 1084883 656740, 108488…
# 3 -0.0857 -0.0292 0.243 0.277 ((1078883 654740, 1080883 654740, 108088…
# 4 -0.0260 -0.0236 0.266 0.292 ((1080883 654740, 1082883 654740, 108288…
# 5 0.212 0.0839 0.331 0.352 ((1082883 654740, 1084883 654740, 108488…
# 6 -0.150 -0.102 0.247 0.265 ((1076883 652740, 1078883 652740, 107888…
# 7 -0.150 -0.102 0.247 0.265 ((1078883 652740, 1080883 652740, 108088…
# 8 -0.0695 -0.0573 0.261 0.280 ((1080883 652740, 1082883 652740, 108288…
# 9 -0.231 -0.194 0.263 0.263 ((1076883 650740, 1078883 650740, 107888…
# 10 -0.190 -0.173 0.280 0.283 ((1078883 650740, 1080883 650740, 108088…
# # ℹ 6,155 more rows
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.
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.
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 .