Geosci. Model Dev., 18, 211–237, 2025
https://doi.org/10.5194/gmd-18-211-2025
© Author(s) 2025. This work is distributed under the Creative Commons Attribution 4.0 License.

HOTSSea v1: a NEMO-based physical Hindcast of the Salish Sea (1980–2018) supporting ecosystem model development
Greig Oldford1,2, Tereza Jarníková, Villy Christensen, and Michael Dunphy
1		Institute for Oceans and Fisheries, University of British Columbia, Vancouver, V6T 1Z4, Canada
2		Ecosystem Sciences Division, Fisheries and Oceans Canada, Nanaimo, V9T 6N7, Canada
3		Tyndall Centre for Climate Change Research, School of Environmental Sciences, University of East Anglia, Norwich, UK
4		Ocean Sciences Division, Institute of Ocean Sciences, Fisheries and Oceans Canada, Sidney, V8L 5T5, Canada
Correspondence: Greig Oldford (greig.oldford@dfo-mpo.gc.ca)
Received: 27 March 2024 – Discussion started: 22 May 2024
Revised: 10 November 2024 – Accepted: 14 November 2024 – Published: 20 January 2025

Published by Copernicus Publications on behalf of the European Geosciences Union.

Abstract. Decadal-scale oceanographic, environmental, and ecological changes have been reported in the Salish Sea, an ecologically productive inland sea in the northeast Pacific that supports the economies and cultures of millions of people. However, there are substantial data gaps related to physical water properties that make it difficult to evaluate trends and the pathways of effects between physical ocean water properties and the productivity of marine ecosystems. With the aim of addressing these gaps, we present the Hindcast of the Salish Sea (HOTSSea) v1, a 3D physical oceanographic model developed using the Nucleus for European Modelling of the Ocean (NEMO) ocean engine, with temporal coverage from 1980–2018. We used an experimental approach to incrementally assess sensitivity to atmospheric and ocean reanalysis products used for boundary forcings and to the horizontal discretisation of the model grid (∼ 1.5 km). Biases inherited from forcings were quantified, and a simple temperature bias correction factor applied at one ocean boundary was found to substantially improve model skill. Evaluation of salinity and temperature indicates performance is best in the Strait of Georgia. Relatively large biases occur in near-surface waters, especially in subdomains with topography narrower than the model grid’s horizontal resolution. However, we demonstrated that the model simulates temperature anomalies and a secular warming trend over the entire water column in general agreement with observations. HOTSSea v1 provided a first look at spatially and temporally heterogenous ocean temperature trends throughout the northern and central part of the domain where observations are sparse. Overall, despite the biases inherited from forcings and a relatively coarse horizontal discretisation, HOTSSea v1 performs well at representing temperature and salinity at the spatial–temporal scales needed to support research related to decadal-scale climate effects on marine ecosystems, fish, and fisheries. We conclude by underscoring the need to further extend the hindcast to capture a regime shift that occurred in the 1970s.
Copyright statement. The works published in this journal are distributed under the Creative Commons Attribution 4.0 License. This license does not affect the Crown copyright work, which is re-usable under the Open Government Licence (OGL). The Creative Commons Attribution 4.0 License and the OGL are interoperable and do not conflict with, reduce or limit each other.
© Crown copyright 2025
shift, within the Salish Sea there are seasonal, annual, and decadal-scale changes to physical oceanographic and atmospheric patterns that have been reported, including changes in seasonal wind patterns (Collins et al., 2009; Masson and Cummins, 2007; Preikshot, 2007; Tuller, 2004), precipitation patterns (Beamish, 1993; Morrison et al., 2002; Yin et al., 1997), ocean water temperatures (Beamish et al., 2010; Masson and Cummins, 2007), properties related to ocean acidification (Feely et al., 2009; Ianson et al., 2016; Jarníková et al., 2022), and river discharge and temperatures (Islam et al., 2019; Martins et al., 2011; Riche et al., 2014). Increasing seasonal stratification and warmer surface waters may have increased the frequency and duration of harmful algal blooms (Esenkulova et al., 2021; Moore et al., 2015). Changes to regional climate patterns appear to have increased the variability of the date of the spring phytoplankton bloom (Allen and Wolfe, 2013), which may have led to spatial–temporal mismatches between predators and prey (Allen and Wolfe, 2013; Suchy et al., 2022) and affected the composition of larval fish assemblages (Guan, 2015). Changing ocean conditions in the Salish Sea are also hypothesised to affect the abundance, composition, and spatial–temporal availability of prey for Pacific salmon via various pathways of effects (Pearsall et al., 2021). Several correlative studies link sea surface temperature and stratification with declining survival of several salmon species in the Salish Sea, particularly juvenile coho salmon (Oncorhynchus kisutch), Chinook salmon (O. tshawytscha), and steelhead (O. mykiss; Beamish, 1995; Pearsall et al., 2021; Perry, 2021; Sharma et al., 2013; Sobocinski et al., 2020, 2021; Walters and Christensen, 2019).
The patchy nature of oceanic data, particularly as we delve deeper into historical records within the Salish Sea, leads to uncertainty about the pace and spatial–temporal patterns of oceanographic change. Sparse observations also limit our ability to detect associations and evaluate mechanistic links between physical oceanographic changes and marine ecosystem dynamics. Physical oceanographic hindcasts are a pivotal tool for addressing such data gaps, offering a retrospective lens through which past oceanic conditions are reconstructed. Physical ocean models may also be coupled or linked to biogeochemical and ecosystem models using an end-to-end approach useful for evaluating mechanistic drivers and dynamic pathways of effects between water properties, marine ecosystems, fisheries, and other human uses (Libralato and Solidoro, 2009; Macias et al., 2014; Piroddi et al., 2021; Rose, 2012; Rose et al., 2010).
Although several oceanographic and biogeochemical models have been developed for the Salish Sea, attributes of these models presently limit their suitability for a long hindcast, including computational cost due to high resolution and a focus on shorter-term simulations (Jarníková et al., 2022; Khangaonkar et al., 2019; Olson et al., 2020; Soontiens et al., 2016; Soontiens and Allen, 2017), too coarse a resolution for use in the Salish Sea due to a focus on the wider BC coast (Peña et al., 2016), or a particular focus on Puget Sound (Khangaonkar et al., 2012, 2019, 2021a; MacCready et al., 2021; Moore et al., 2015). It is therefore one aim of this study to develop a physical hindcast with a particular focus on the central and northern portion of the Salish Sea (i.e. the Strait of Georgia) with adequate spatial–temporal resolution to enable a long hindcast and acceptable model skill for supporting marine ecological research and ecosystem management. Although the acceptable model error will depend on the specific research question, ecological patterns and associated processes related to plankton and fish are often orders of magnitude greater than those required to study physical and chemical processes in marine ecosystems (Fulton et al., 2019). As a local example, hindcasting the timing of the spring phytoplankton bloom in the Strait of Georgia within an error margin of several days to 1 week at subregional scales is useful for studying spatial–temporal mismatches in predators and prey (Allen and Wolfe, 2013; Gower et al., 2013; Suchy et al., 2022). Many factors hypothesised to be mechanistically linked to growth and mortality of juvenile salmon exhibit variability on timescales of 1 week or more and at spatial scales of >10 km (Pearsall et al., 2021). Another hinderance to development of a long hindcast for the area has been a lack of atmospheric and oceanic data extending back to at least 1980 to use as model forcings, and there are doubts about the adequacy of the few existing products with respect to spatial–temporal resolution. A second aim of this study is therefore to conduct experimental evaluation to identify biases inherited from external forcings (including some that only recently have been made available) to determine to what degree the available forcings limit the development of long oceanographic hindcasts in the area.
Here, we present HOTSSea v1, developed using the Nucleus for European Modelling of the Ocean (NEMO) Ocean engine (Madec et al., 2017). We describe and give a rationale for the model setup with attention to three aspects that are particularly important for developing a long hindcast for the domain: (1) the biases inherited by using various atmospheric and ocean reanalysis products as surface and boundary forcing; (2) the effect of applying temperature bias corrections to the open-ocean boundary forcing; and (3) a preliminary assessment of model performance relevant to the aforementioned research applications, including decadal-scale trends. Priority areas for improvement and further evaluation are also highlighted, and, finally, we use the model to provide a first look at decadal-scale trends in the central and northern portion of the domain where historical observations are especially sparse.
2		Model overview
The NEMO ocean engine, version 3.6, supports simulations of ocean dynamics and thermodynamic processes in three dimensions (Madec et al., 2017). The physical model framework is governed by primitive equations under hydrostatic balance using the Boussinesq approximation where density variations are neglected except in their contribution to the buoyancy force (Madec et al., 2017). HOTSSea v1 was implemented in a high-performance computing cluster (Digital Research Alliance of Canada, 2022), and the model’s scope is limited to physics and hydrodynamics (e.g. tides, salinity, temperature) – biogeochemistry is not included in HOTSSea v1 due to computational cost. We used state variables of practical salinity (PSU) and potential temperature (°C), as well as the EOS-80 equation of state (Millero, 2010). Sea ice is not included in HOTSSea v1 given that it occurs only occasionally in deep inland waters of fjords such as Jervis Inlet. To address issues of omitting ice, we applied NEMO’s ice-if option, where the water temperatures are limited to the local salinity-dependent freezing point.
2.1		Spatial–temporal configuration
The model domain of HOTSSea v1 includes several distinct geographic areas within the Salish Sea: the Juan de Fuca Strait, Strait of Georgia, Gulf Islands, and Puget Sound (Fig. 1). A key application of the model will be to provide forcings for biogeochemical and ecosystem models developed to investigate decadal-scale change. The model’s spatial domain therefore fully encompasses the domain of an ecosystem model under parallel development that focuses on the Strait of Georgia and uses the Ecospace model framework (de Mutsert et al., 2023; Walters et al., 1999). To ensure the dynamics of the Strait of Georgia were resolved, it was deemed important to resolve the broader Salish Sea, with the connection to the open ocean via the Gulf Islands and the Juan de Fuca Strait being particularly important for resolving tides and estuarine flow (Ebbesmeyer et al., 1989; MacCready et al., 2021). The horizontal grid used in NEMO is discretised on a curvilinear orthogonal Arakawa C grid generalised to three dimensions (Arakawa and Lamb, 1977; Madec et al., 2017). The basic spatial–temporal configuration of HOTSSea v1 began with a previous configuration, SalishSeaCast, implemented at approximately 500m horizontal resolution for the same domain (Olson et al., 2020; Soontiens et al., 2016; Soontiens and Allen, 2017). The ∼ 500 m horizontal resolution grid and bathymetry used in the SalishSeaCast model were reduced by a factor of 3 in each horizontal direction, taking the mean depth of the neighbouring cells to assign the new depths. The new grid is approximately 1.5 km in horizontal resolution and has a width of ∼ 200 km and length of ∼ 450 km (132 cells × 299 cells; Fig. 1), which matches the grid and resolution of the Ecospace ecosystem model. The grid is rotated 29° anticlockwise to true north to align with the axis of the Strait of Georgia. The bathymetry was already processed for Salish- SeaCast to avoid sudden changes in depths across grid cells and maintain open channels in narrow passages. We made additional manual edits to maintain channels between islands,

Figure 1. Map of model domain showing geographic features, extents of the HOTSSea NEMO model domain (medium grey), and bathymetry.

maintain connectivity of the main Fraser River channel to the outflow, and avoid erroneously isolating bodies of water. Some narrow water bodies such as Sechelt Inlet, Salmon Inlet, Burrard Inlet, and the Indian Arm fjord are not resolved in this setup (outlines of these areas are visible in Fig. 1). The depths of edited channel cells were approximated from depth averages taken from ∼ 80m resolution bathymetric data (PSF, 2014). To ensure tidally driven dynamics were not lost, the main channel of the Fraser River was extended inland by manually adding non-existent river channel cells approximately 150 km in total length, following Soontiens and Allen (2017).
The vertical grid for HOTSSea v1 is divided into 40 vertical (z) levels that are gradually stretched to achieve higher resolution at the surface, ranging from 1 m vertical resolution in the upper 10 m to approximately 27 m widths at the deepest level (420 m). Partial steps were enabled to limit large changes in bathymetry between adjacent grid cells. The thickness of each layer is proportionally scaled at each time step as sea surface height (SSH) changes using a non-linear free-surface scheme referred to as the “variable volume option” (Levier et al., 2007). HOTSSea v1 uses a non-linear free-surface option to time-split the solving of the barotropic and baroclinic free surface. The barotropic and baroclinic time steps are set to 6 and 120 s, respectively, and the vertical momentum and tracer advection time stepping is set to 2 s. The model was run from 1 January 1979 to 1 January 2019, where the 1980 atmospheric forcings were duplicated and applied to 1979, such that we treat 1979 as a model spin-up year and exclude it from evaluation. A 1-year spin-up was based on a minimum estimate of deepwater residency time which elsewhere has been reported to range between 1 and 3 years (Pawlowicz et al., 2019). Initial conditions for January 1979 for temperature and salinity across the domain were generated using climatologies for December and January using SalishSeaCast outputs from 2007 to 2020. An experimental bias correction to the Ocean and Sea Ice Reanalysis v5 (ORAS5) temperature fields was applied when running the final hindcast.
2.2		Boundary conditions and forcings
2.2.1		Atmospheric
The Regional Deterministic Reforecast System (RDRS v2.1; Gasset et al., 2021) supplied the atmospheric conditions for forcing the full HOTSSea v1 hindcast. RDRS v2.1 is currently the highest-resolution atmospheric reanalysis product available extending back to 1980 (0.09°; ∼ 10 km horizontal). Two additional atmospheric forcings (Table 1) were evaluated as part of an experimental design: the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5, a global reanalysis product extending back to 1979 (hourly, at approximately 31 km horizontal resolution; Dee et al., 2011; Hersbach et al., 2020), and the High Resolution Deterministic Prediction System (HRDPS), with spatial coverage of the northern part of North America (Canada and northern United States) and with hourly coverage at ∼ 2.5 km horizontal resolution for 2014–2020 (ECCC, 2020). The RDRS v2.1 product occupies an intermediate horizontal resolution between ERA5 and HRDPS, and the three together offered an opportunity to explore the effect of horizontal resolution of atmospheric forcing on model performance in the Salish Sea.
2.2.2		Open boundaries
There are two boundaries that connect the Salish Sea to the Pacific Ocean: the mouth of the Juan de Fuca Strait in the southwest and Johnstone Strait in the north (Fig. 1). To first evaluate the effects of using different ocean boundary forcings at the mouth of the Juan de Fuca Strait, a higher- resolution model, CIOPS-W (Lin et al., 2022), was used in shorter evaluation runs (horizontal resolution 2–2.5 km; 1/36°; Tables 1 and 2). ORAS5 (Tietsche et al., 2017; Zuo et al., 2019) was the only available reanalysis product with coverage for the full model hindcast and was used to supply open-ocean boundary conditions in the final model. ORAS5 has a horizontal resolution at the latitude of the Salish Sea of approximately 18 km (0.25°). At the northern boundary (Johnstone Strait), we used a monthly climatology of temperature and salinity (Dosser et al., 2020, 2021).
2.2.3		River discharge and runoff
River input into the Salish Sea periodically creates a brackish layer extending across the Strait of Georgia and drives strong estuarine circulation via the Juan de Fuca Strait (Harrison et al., 1983). The Fraser River is the largest single source of freshwater influx into the domain and supplies approximately two-thirds of the total annual freshwater input (Pawlowicz et al., 2019). Fraser River discharge is monitored as part of a long-term programme (Morrison et al., 2012). Following Soontiens and Allen (2017), we used available flow records for the Fraser River from gauges approximately 150 km inland at the city of Hope, BC (DFO, 2024), and supplemented the Fraser River flow data with climatological data for additional freshwater input downstream of the station. A climatology was used for Fraser River runoff temperatures (Morrison et al., 2002) due to a lack of long-term measurements from the lower Fraser River. The location of river outflow for the Fraser River was placed in the main channel before the river branches into a delta (at the town of Delta, BC). All other river outflows were assigned to the grid cell closest to the river mouth. Many rivers other than the Fraser River are not monitored, so climatological patterns for discharge and temperature for 150 rivers flowing into the Salish Sea were used (Morrison et al., 2012). We adapted the input file containing these river input data from the ∼ 500 m horizontal resolution model grid used by Soontiens et al. (2016) to the ∼ 1.5 km horizontal resolution used here and adjusted the outflow locations as required.
2.2.4		Tides
At the two open boundaries, tides were forced with eight tidal constituents (K1, O1, P1, Q1, M2, K2, N2, and S2). Tidal heights and currents at the Juan de Fuca Strait boundary were originally taken from WebTide (Foreman et al., 2000) and then manually tuned (Soontiens et al., 2016). At the northern open boundary sea surface height and tidal harmonics were forced for the major M2 and K1 constituents, and SSH harmonics for the O1 and S2 harmonics were configured using calculations from Thomson and Huggett (1980), with remaining constituents taken from WebTide and subsequently tuned.
3		Model evaluation
Observations were collated from various instruments and sources (Table 2) and used to do a preliminary evaluation of the model’s performance with respect to sea surface temperature (SST), sea surface salinity (SSS), and temperature and salinity over depths. To understand the trade-offs between spatial–temporal resolution, tractability, and model
Table 1. External forcing used in the model (n/a: not applicable).
Forcing dataset	Forcing type	Model runs	Temporal extent and resolution	Horizontal resolution	Citation	
Regional Deterministic Reforecast System (RDRS v2.1)	Surface/atmospheric	Final	1980–2018; hourly	0.09°; ∼ 10km	Gasset et al. (2021)	
European Centre for Medium-Range Weather
Forecasts (ECMWF) Ocean and Sea Ice Reanalysis v5 (ORAS5)	Open-ocean boundary conditions	Final	1975–2018;
monthly	0.25°; ∼ 18km	Tietsche et al. (2017);
Zuo et al. (2019);
Copernicus Climate
Change Service (2021)	
Runoff/river climatology and gauge data	Runoff	Final	1979–2018; hourly and daily	n/a	Morrison et al. (2012);
Soontiens et al. (2016);
Water Survey of Canada (2023)	
Tidal constituents	Tidal forcing at open boundaries	Final	n/a	n/a	Soontiens et al. (2016)	
Coastal Ice Ocean Prediction System (CIOPS) West	Open-ocean boundary conditions	Evaluation	2007–2019; hourly	1/36°; ∼ 2.5 km	Lin et al. (2022)	
ECMWF ERA v5 (ERA5)	Surface/atmospheric	Evaluation	1979–present; hourly	0.28°; ∼ 31 km	Dee et al. (2011); Hersbach et al. (2020)	
High Resolution	Surface/atmospheric	Evaluation	2014–2020; hourly	0.0225°; ∼ 2.5 km	ECCC (2020)	
Deterministic Prediction
System (HRDPS)

skill, we used an experimental approach where forcings were incrementally swapped to help with isolating the most likely source of model error and bias. The NEMO-based Salish- SeaCast model (v201905) outputs were used for comparison when evaluating the effect on overall model performance of changing the spatial–temporal setup, as we used this model as a foundation for the HOTSSea v1 model. As such, HOTSSea v1 shares limitations of the SalishSeaCast model, such as no wetting-and-drying capability in intertidal areas, climatologies used for river flow for all rivers except for the Fraser River, and biases in temperature and salinity. Aspects of SalishSeaCast’s model skill with respect to physical properties have been previously reported (Olson et al., 2020; Soontiens et al., 2016; Soontiens and Allen, 2017). The mean temperature biases over all depths ranged from -0.21 to 0.13 °C, and for depths less than 15 m they ranged from -0.34 to 0.36 °C, respectively (Olson et al., 2020; Soontiens et al., 2016; Soontiens and Allen, 2017). The mean salinity bias over the same two depth strata ranged from -0.74 to +0.23 PSU and -1.62 to 0.23 PSU, respectively (Olson et al., 2020; Soontiens et al., 2016; Soontiens and Allen, 2017). Our results may differ because the previous study did not necessarily use the same observational data nor did it use the subdomain definitions we used. For evaluation of temperature trends in the final hindcast, we used the modelled longterm temperature trend against observations at Nanoose station, the only long-term dataset with at least biweekly depth profiles done in the model domain extending back to the beginning of the hindcast (Table 2).
3.1		Experimental evaluation
The years chosen for running preliminary experimental evaluations were 2016–2018. The available forcing data and models had coverage for those years, and generally a larger volume of evaluation data is available for the most recent years (Table 3). Experimental runs of the HOTSSea model were given run codes. The first was HOTSSea v0.1, which used the highest-resolution atmospheric and ocean boundary forcings available. This run was used to do a comparative evaluation with SalishSeaCast v201905, which has a higher horizontal resolution (∼ 500 m versus ∼ 1500 m). The LiveOcean model (Fatland et al., 2016) forcings used in SalishSeaCast at the Juan de Fuca Strait (JFS) open-ocean boundary were not available for the 2016–2018 period, so in HOTSSea v0.1 we used CIOPS-W BC12, a model also developed using the NEMO v3.6 ocean engine covering the northeast Pacific at an approximate horizontal resolution of 2.0– 2.5 km (Lin et al., 2022). The HOTSSea v0.12 experiment was used to evaluate the effect of swapping from HRDPS to the ERA5 atmospheric forcings (∼ 31 km horizontal; Dee et al., 2011; Hersbach et al., 2020). At the time HOTSSea development began in 2021, ERA5 was the only climate reanalysis product available for the entirety of the hindcast period. At the time of writing, it is still the only reanalysis extending back to the 1940s, and therefore evaluation of this product
Table 2. Summary of data used for model evaluation.
Instrument type	Dataset title	Variables	Observations
(N)	Description	Source	
Conductivity, temperature, and depth (CTD) casts	Fisheries and Oceans Canada’s (DFO) Institute of Ocean Sciences (IOS) CTD casts dataset	Conductivity, temperature, depth, pressure, oxygen and salinity	24 810	Contains CTD measurements collected in the central Strait of Georgia, British Columbia, Canada, using rosette-mounted CTDs.	DFO (2022b)	
	DFO IOS (including Nanoose station)	Salinity, temperature, depth, pressure	3942	Surveys conducted from 1965 to present including Nanoose Bay station, a Canadian military CTD dataset, which were provided upon request from DFO.	https://waterproperties.
ca (IOS, 2025)	
	Hakai Institute	Salinity, temperature, depth, pressure	2871	CTD data collected from 2012 to present by the Hakai Institute in waters surrounding Calvert Island, Johnstone Strait, and Quadra Island areas.	Hakai Institute (2024)	
	Pacific Salmon Foundation (PSF)	Salinity, temperature, depth, pressure	3437	CTD casts collected by PSF for the Strait of Georgia.	PSF (2023)	
Lightstation (LS) nearsurface water properties	Environment and Climate Change Canada (ECCC)	Temperature, salinity	7	Observations from lightstations where daily sea surface temperature and salinity measurements have been collected from 1914 to present. Measurements were made daily using seawater collected in a bucket lowered into the surface water at or near the daytime high tide.	DFO (2022a)	
Buoys	ECCC via DFO	Sea surface temperature (SST)	5	Wave and temperature data from buoys. Sea surface temperature data have undergone automated quality control. Historical data are merged with real-time acquisition.	Fisheries and Oceans
Canada (DFO) (2024)	

for use for atmospheric forcing was a priority. The v0.14 and v0.16 experiments jointly helped evaluate the effect of using the ORAS5 (∼ 18 km horizontal; Tietsche et al., 2017; Zuo et al., 2019) dataset for ocean boundary conditions at the mouth of the Juan de Fuca Strait, which will be used in the final hindcast; the two experiments used the lowest-resolution atmospheric forcings (v0.14; ERA5) and the highest-resolution atmospheric forcings (v0.16; HRDPS) to assist with isolating the effects on model performance of ORAS5 versus the atmospheric forcings. The HOTSSea v0.18 experiment used the RDRS v2.1 atmospheric outputs for forcing, which have an intermediate horizontal resolution of ∼ 10 km (Gasset et al., 2021). To evaluate each experiment, we used data, methods, and statistics as described in the next section. Model performance was evaluated using results aggregated over the 2016–2018 period – analyses were also carried out on model results grouped by month and year, though only results aggregated for the entire period are presented here, and only the results using CTD measurements are highlighted here for brevity.
3.2		Model–observation evaluation methods
The final hindcast was evaluated using the datasets grouped by subdomains (Table 2, Fig. 2). Subdomains were selected based on distinct geographic features, data availability, and physical characteristics.
3.2.1		Vertical profiles
CTD casts were acquired from various sources (Table 2). After quality control, 27 272 CTD casts were used in the analysis of the final hindcast. These data have heterogeneous spatial coverage when aggregated by subdomain (Fig. 2; Table 4). For the model intercomparison and experimental evaluation only CTD data from 2016–2018 were used (N = 8012), which also had spatially heterogenous coverage across subdomains.
The closest model grid cell and time index was found for each CTD measurement, and the measurements over depths for each CTD cast were vertically interpolated to the model depth levels. To calculate the statistics described below, CTD data were first grouped by subdomain, period, and depth strata during analysis. We highlight in the subsequent sections only the results of grouping the data first by subdomain and model depth level and second by subdomain and selected depth strata (0 to >30 m, 30 to >150 m, >150 m, and over all depths). For statistics grouped by depth strata, the depth-integrated mean from individual CTD casts within each depth grouping was first calculated. These values were treated as a single measurement (oi ) in the set of CTD casts,

Table 3. Experimental evaluation run codes and forcings.
Model	Version or run code	Evaluation purpose	Years	Surface forcing	Ocean boundary forcing	Reference	
SalishSeaCast (∼ 0.5 km)	v201905	Comparison of model performance at higher horizontal resolution	2016–2018	HRDPS	LiveOcean at JFS boundary, climatology at northern boundary	Soontiens et al. (2016);
Soontiens and
Allen (2017);
Olson et al. (2020)	
CIOPS-W	BC12	Comparison of model performance at lower horizontal resolution	2016–2018	HRDPS and
RDPSa (∼ 10 km)b	Regional Ice Ocean
Prediction System
(RIOPS) v2	Lin et al. (2022)	
HOTSSea (∼ 1.5 km)	v0.1	Comparison with two models listed above using ∼ 1.5 km2 horizontal resolution	2016–2018	HRDPS	CIOPS-W	This study	
HOTSSea	v0.12	Evaluates sensitivity to lower-resolution atmospheric forcing	2016–2018	ERA5	CIOPS-W	This study	
HOTSSea	v0.14	Evaluates sensitivity to lower-resolution ocean boundary and atmospheric forcings	2016–2018	ERA5	ORAS5 (western boundary) and climatology (northern boundary)	This study	
HOTSSea	v0.16	Uses the highest-resolution atmospheric forcings available; evaluates sensitivity to lower-resolution ocean boundary forcings	2016–2018	HRDPS	ORAS5 (western boundary) and climatology (northern boundary)	This study	
HOTSSea	v0.18	Uses the forcings used in the
final model; evaluates sensitivity to intermediate-resolution atmospheric
and lower-resolution boundary forcings	2016–2018	RDRS v2.1	ORAS5 (western boundary) and climatology (northern boundary)	This study	
HOTSSea	v1.01	First full hindcast run	1979–2018	RDRS v2.1	ORAS5 (western boundary) and climatology (northern boundary)	This study	
HOTSSea	v1.02	Full hindcast run with bias correction	1979–2018	RDRS v2.1	ORAS5 – temperature bias correction (western boundary) and climatology (northern boundary)	This study	
a RDPS: Regional Deterministic Prediction System (used by CIOPS-W; not used herein). b Approximate horizontal resolutions are included in brackets where not otherwise included in Table 1.

Table 4. Spatial distribution of CTD data used for evaluation of the entire hindcast (1980–2018) and shorter experiments (2016–2018).
Subdomain	CTD count∗	
	1980–2018	2016–2018	
Discovery Islands (DI)	3649 (2.3 km-2)	1884 (1.2 km-2)	
Strait of Georgia North (SGN)	12 365 (2.9 km-2)	3080 (0.7 km-2)	
Strait of Georgia South (SGS)	4140 (1.5 km-2)	1517 (0.5 km-2)	
Gulf Islands (GI)	2512 (1.4km-2)	871 (0.5 km-2)	
Haro Strait (HS)	800 (2.0 km-2)	184 (0.5 km-2)	
Puget Sound (PS)	99 (0.04 km-2)	23 (0.01 km-2)	
Juan de Fuca Strait (JFS)	3707 (0.9 km-2)	453 (0.1 km-2)	
Total:	27 272 (1.5 km-2)	8012 (0.4 km-2)	
∗ Brackets indicate spatial density of measurements.

Figure 2. Map of HOTSSea v1 model domain (red rectangle) and subdomains used for analysis (black polygons; DI, Discovery Islands; SGN, Strait of Georgia North; SGS, Strait of Georgia South; HS, Haro Strait; GI, Gulf Islands; JFS, Juan de Fuca Strait; PS, Puget Sound). Locations of CTD casts are indicated by orange stippling (darker denotes higher density), and Nanoose station is indicated by the black star.

Lightstation ID	Years	Location (subdomain)	
Active Pass	1980–2011	SGS	
Cape Mudge	1980–1985	SGN	
Chrome Island	1980–2019	SGN	
Entrance Island	1980–2018	SGN	
Race Rocks	1980–2018	JFS	
Sheringham Point	1980–1988	JFS	
Sisters Islets	1980–2008	SGN	

Table 5. Lightstation data summary used for evaluation of SST and SSS.

(2)

N , across each subdomain, with o representing the mean of depth-integrated means. Model results were extracted for each observation and depth-integrated in the same manner.
The model error was calculated for each model– observation pair (m, o) of time series (error = m - o) with the bias being the mean error. The root mean square error (RMSE) was also calculated for each depth stratum, time frame, and subdomain. The bias and RMSE are often used to infer the accuracy of a model, whereas the centred root mean square error (CRMSE) quantifies the precision (Walther and Moore, 2005) as the variability of the model as compared with observations as
u1N
CRMSE =	((mi-m) - (oi-o))2,	(1)
N i=1
where oi denotes a single observation out of a set of observations, N, with the corresponding model output denoted by mi. The Willmott skill score (WSS) was also used, a dimensionless measure of model skill ranging from zero, denoting poor agreement between model and observations, to one, denoting perfect agreement (Willmott, 1981):
PN i=1(mi-oi)2
WSS = 1 N	.
iN=1(|mi -o| + |oi - o|)2
3.2.2		Sea surface temperature and salinity
Sea surface temperature (SST) and salinity (SSS) were evaluated using measurements collected at high tide during daylight hours by lighthouse staff at seven lighthouses throughout the domain (Fig. 2; Treasury Board Secretariat, 2023). Lightstation data (LS) were available for the entirety of the hindcast period for Chrome Island, Entrance Island, and Race Rocks lighthouses, whereas others had partial coverage (Table 5). The time of day when samples were taken was not always provided in the dataset. As such, the closest tidal gauge each day was found for each lighthouse, and the high-tide time was extracted and then matched to the lighthouse sample data as required.
Sea surface temperature measurements were also available from buoys within the model domain (Fig. 2). Canadian buoy

R=

11 σoσm N

(3)

(4)

(5)

(6)

Figure 3. Monthly climatology of ORAS5 temperature bias over depths at the Juan de Fuca Strait open-ocean boundary.

Table 6. Buoy data summary used for evaluation of SST.
Buoy ID	Years	Location (subdomain)	
46131	1992–2019	SGN	
46134	2001–2019	GI	
46146	1992–2019	SGN	

data were downloaded from online repositories (Fisheries and Oceans Canada (DFO), 2024; NOAA National Buoy Data Centre, 1971). The buoy data we prepared had heterogeneous temporal coverage of the hindcast period (Table 6). For the evaluation of the full hindcast we only used buoys with 10 or more years of data (buoy IDs: 46146, 46131, and 46134).
The statistics used for evaluation of SST and SSS include the standard deviation, σ , and Pearson’s R (i.e. the correlation coefficient):
N
(mi - m)(oi - o).
i=1
For plotting of results using Taylor diagrams, the σm and CRMSE statistics were normalised by scaling to the standard deviation of the observations to allow the display of evaluation results for multiple stations on the same plots. The normalised centred root mean square error (NCRMSE) is computed using the relationship between σ , R, and CRMSE,
CRMSE = σo2 + σm2 -2σoσmR;
thus,
NCRMSE = CRMSE / σo,
and
0	σm σm = m	σo where σm0 denotes the normalised standard deviation. Note that in target plots, the NCRMSE has been modified by the sign of σm - σo (Kärnä et al., 2021).
3.3		Ocean boundary temperature bias correction
The ORAS5 global ocean reanalysis extends back to 1958 (Copernicus Climate Change Service, 2021), and it would be ideal to extend the model back to 1958, too; however, the magnitude of biases that would be inherited from ORAS5 was hitherto unknown. Based on the results of experiments (Table 3), we suspected biases could be substantial. To investigate, we collated observations from CTD instruments sampled between 1980 and 2018 from the mouth of the Juan de Fuca Strait (N = 2162) and compared these observations to temperatures in the ORAS5 data. We chose to do this for
the entire hindcast period rather than using the experiments for only 2016–2018 given the observational data available at this location for the experimental period were relatively limited. ORAS5 outputs were interpolated from monthly to daily using the cdo toolset (Schulzweida, 2022), which is the same procedure done by NEMO internal routines during model runs. Each CTD measurement was then matched to the closest ORAS5 grid point. Both the observations and the ORAS5 model data were interpolated vertically to the HOTSSea model depth levels. The ORAS5 model bias was calculated using monthly mean bias at each depth level.
The analysis indicated that ORAS5 at the JFS boundary is biased warm at most depths and months except January, with the mean summer bias near the surface approaching +4 °C (Fig. 3). Surface waters in the ORAS5 model were biased fresh in the spring and summer, especially in the top 10m where the mean bias approached -4 PSU. As a first step towards a more comprehensive bias correction, the mean monthly temperature bias for each depth level was used as a correction factor applied to the boundary conditions. We chose to prioritise temperature to isolate a single variable and because applying a bias correction factor to salinity would run a risk of introducing dynamic instability. The temperature bias correction method is acknowledged here to be a crude approach compared to various alternatives (Adachi and Tomita, 2020). The model run with the ocean boundary temperature bias correction factor applied to ORAS5 is referred to as HOTSSea v1.02. We then evaluated the change in model performance versus with no temperature bias correction (HOTSSea v1.01).
3.4		Trend analysis
The ability of the HOTSSea model to recreate observed longterm ocean temperature trends is one of the most important tests of the utility of the model. The only station in the model domain where measurements over the entire water column were collected with biweekly regularity for the entirety of the hindcast is Nanoose station, located in the central Strait of Georgia (Fig. 2). At this location, CTD casts have been sampled approximately every 2 weeks since 1979 and with
less regularity back to the late 1960s. Trends in water temperature at Nanoose station were previously analysed and estimated to be 0.24 ± 0.1 °C per decade between 1970 and 2005 (Masson and Cummins, 2007; MC07). We acquired the same Nanoose station data as analysed by MC07, with updates for recent years from a digital archive at the DFO Institute of Ocean Sciences (IOS) and in the DFO IOS Water Properties Database in August 2022 (IOS, 2025). To cross-check our analysis with MC07, we also calculated the slope for the same period as that study (1970–2005), similarly using ordinary least squares linear regression, and found a similar result: a warming trend over the water column (4.5–400 m) of 0.257 °C per decade was assessed here versus 0.24 °C per decade in MC07. We attribute the difference compared to the application of a boxcar filter to the anomalies in the previous study.
To compare temperature anomalies and trends observed at Nanoose station with those predicted by the model, first obviously erroneous measurements (e.g. with coordinates on land or depths exceeding the maximum depth of ∼ 400 m at Nanoose station) were removed and units were converted to match those used in HOTSSea. The number of CTD casts ultimately used was 5692. Measurements taken in each CTD cast were interpolated to match the model depth levels of the closest HOTSSea grid cell. The mean for each depth level grouped in 2-week time intervals was calculated for the hindcast period. A climatology over depth and time was generated for the hindcast period by taking the mean across years for each depth level and 2-week block. Time series of temperature and salinity anomalies were calculated by differencing the depth-binned data and the climatology. A blind approach was then used to match each CTD cast to HOTSSea outputs where model data were only extracted for dates, times, and locations corresponding to each CTD measurement (i.e. gaps in the Nanoose time series were also present in the extracted model results). The same procedure as above was then used for the model results to generate climatologies and anomalies.

TS = median

xi - xj
j-i

for all i < j ,

(7)

Analysis of temperature trends was conducted using the Nanoose station data and model results drawn from the same area to examine the performance of HOTSSea v1.0x at simulating long-term trends. Depths shallower than 4.5 m were omitted because the first depth at which measurements were taken was inconsistent (following MC07). The analysis of the depth-integrated anomaly trend over the entire water column (4.5–400 m) was first done, followed by analysis of the trend at each model depth. To quantify the magnitude of a linear trend, the Theil–Sen slope was used (TS; Sen, 1968; Theil, 1950). The TS approach is more robust to data gaps, outliers, and non-normal and heteroscedastic residuals than ordinary least squares linear regression (LR; Wilcox, 1998). To calculate TS, the median slope of all data pairs is found: where (xi , xj) is a pair of values in the ordered time series (j >i). The LR method used in MC07 and the TS method were compared, and we estimated the trend for the 1970– 2005 period using TS to be 0.256 °C per decade, effectively the same as the LR method after accounting for measurement precision (see above). The similarity between trends estimated between the LR and TS method is consistent with findings in a previous analysis of SST in the same region (Amos et al., 2015).
Detrended residuals were analysed for periodicity, autocorrelation, non-normal distribution, and heteroskedastic- ity prior to choosing the test used for determining statistical significance. The data were first “deseasoned” using the climatological mean to produce anomalies, as described above. The deseasoned anomalies were detrended by removing the secular trend calculated using the TS method, and fast Fourier transform (FFT) was used to detect any remaining periodicities in the detrended and deseasoned residuals by examining the top five peaks in the power spectrum (Bluestein, 1970; Cooley and Tukey, 1965) using scipy (Virtanen et al., 2020) and numpy Python packages (Harris et al., 2020). The presence of autocorrelation was evaluated by modelling the residuals as a first-order autoregressive process (AR-1) and computing the autocorrelation function (ACF) using Bartlett’s method for computing the 95 % confidence interval (Brockwell and Davis, 2016; Parzen, 1964) with the statsmodels package in Python (Seabold and Perktold, 2010). The residuals were tested for normality using the Shapiro– Wilk test, where a p value of ≤ 0.9 was interpreted as insufficient evidence to reject a non-normal distribution (Shapiro and Wilk, 1965) using scipy. Heteroskedasticity was evaluated using White’s test (White, 1980) and the Goldfeld– Quandt test (Goldfeld and Quandt, 1965) using statsmodels.
The analysis using FFT revealed the presence of periodicity in residuals of frequencies of 5.1 and 17.7 years, confirming that the data were deseasoned at the sub-annual scale but suggesting climate modes operating at longer timescales may limit the ability of accurate secular trend estimation, a similar result to previous studies (Amos et al., 2015). The 17.7 periodicity detected is consistent with the Pacific Decadal Oscillation (PDO), the North Pacific Gyre Oscillation (NPGO), or a combination of the two. The 5-year period is approximately consistent with the El Niño–Southern Oscillation (ENSO) cycle (Andres Araujo et al., 2013; McPhaden et al., 2006).
Analysis of the data revealed that the residuals of the detrended anomalies (both model and observation) had significant autocorrelation, were not normally distributed, and had some heteroskedasticity detected at some depths. The Mann– Kendall (MK) method (Kendall, 1948; Mann, 1945), a nonparametric method, was chosen for evaluating trend significance. This method is widely used in climatological and meteorological applications (Amos et al., 2015; Gocic and Tra- jkovic, 2013) and has the advantage of being relatively robust to non-normal data distributions and the data gaps present in the Nanoose time series. To check for statistical significance
using MK, first the S statistic is calculated,
n-1 n
S = X X sgn(xj -xi) ,	(8) i=1 j=i+1
where n is the number of data points in the time series and sgn(xj - xi ) represents the sign function,

sgn(θ) =

+1, ifθ > 0;
0, if θ = 0;	. -1, ifθ < 0 

(9)

When the S statistic is greater than 0, it indicates that values earlier in the time series tend to be lower than those later, and the trend therefore tends positive. Next, the variance of the test statistic is required, which accounts for ties in the data:
1
VAR(S) =	[n(n - 1)(2n+5) 18
- X tp (tp - 1)(2tp + 5) , p=1

(10)

where n is the length of the time series, m is the number of tied values, and tp is the number of ties in the pth tied value. Using S and VAR(S), one computes the standardised normal variate Z:

ZMK =

√ S-1	, if S > 0;
VAR(S) ,
0, if S = 0;
S+1
√VAR(S) , if S < 0

(11)

The standardised statistic, ZMK , follows the normal distribution with E(Z) = 0 and V (Z) = 1. The null hypothesis – that there is no trend – may be rejected if the absolute value of Z is larger than the theoretical value of Z1-α/2 for a onetailed test or Z1-α for a two-tailed test (used here), where α is a chosen statistical significance level, which here was 95 %. The Z values were also used for lower confidence limits (LCLs) and upper confidence limits (UCLs). The MK test, similar to LR with a t test, requires the data be independent, and therefore serial autocorrelation must be dealt with in advance (Helsel and Hirsch, 1992). Yue et al. (2002) proposed a “pre-whitening” procedure to remove the effect of serial autocorrelation prior to estimating the trend significance using MK. When autocorrelation was present, we used the 3PW algorithm (Collaud Coen et al., 2020), which combines three pre-whitening methods to minimise risks of type I and type II errors, which we applied using the mannkendall (v1.1.1) Python code.
After comparing the model with the measurements from the Nanoose dataset, we evaluated the trend in temperature and salinity for several depth strata (<30, 30–150, >150 m, and all depths) in each grid cell in the subdomain using the same statistical methods described above, though the 3PW method was not applied due to computational cost. The analysis was limited to the Strait of Georgia and surrounding waters for several reasons: (i) preliminary results indicated that model performance was best in this part of the domain, (ii) collated observations have relatively more coverage in this area, (iii) and it is the focus of an ecosystem model under development in parallel. Mean weekly water temperatures were calculated for each grid cell, season (winter: December, January, and February; spring: March, April, and May; summer: June, July, and August; fall: September, October, and November), and depth group. When generating plots, grid cells that were shallower than a threshold set for each depth strata were masked (thresholds: 20 m for 0–>30 m; 150 m for 30–>150 m; 200 m for >150 m) and coloured grey.
4		Results
4.1		Experimental evaluation
In the first model experiment (years 2016–2018), HOTSSea v0.1 was compared with the higher-resolution SalishSeaCast model. A notable difference was that v0.1 had a mean bias taken over all depths of 0.14 °C in the Strait of Georgia South (SGS) subdomain versus -0.017 °C in SalishSeaCast (Fig. 4c). Large differences were noted in the northernmost Discovery Islands (DI) subdomain where near-surface (0– 1.5 m) biases were larger in v0.1 relative to SalishSeaCast, approaching +2 °C and -4 PSU at 0.5 m (Fig. 4a, d). We attribute the warm bias in the DI subdomain to the 3-times- coarsened HOTSSea model grid relative to SalishSeaCast, limiting the model’s ability to resolve the relatively narrow topography in the DI. At depths >1.5 m HOTSSea v0.1 performed well, even in the DI, with mean biases over the water column of -0.12 °C and -0.4 PSU. In the JFS subdomain, HOTSSea v0.1 performed similarly to SalishSeaCast except for the introduction of a temperature bias of 0.3–0.5 °C at depths >2 m, with the bias most pronounced at depths >20 m (Fig. 4b). A similar bias was also present in the Strait of Georgia South (SGS) subdomain (Fig. 4c). The boundary forcings were hypothesised to be the main culprit, since different boundary forcings were used here (CIOPS-W) than in SalishSeaCast (LiveOcean).
The v0.12 and v0.14 experiments were designed to evaluate the effect on the model’s performance of using relatively coarse atmospheric forcings from ERA5 (∼ 31 km horizontal resolution) versus the relatively highly resolved HRDPS (∼ 2.5 km horizontal) forcings used in v0.1. The v0.12 and v0.14 experiments used different ocean boundary forcings, though the resulting salinity and temperature biases over depths were indistinguishable between the two (which prompted us to remove v0.12 from Fig. 4). The similarities suggest that any performance impacts of using ORAS5 (∼ 18 km horizontal resolution) versus CIOPS-W (∼ 2.5 km horizontal resolution) were masked by the relatively greater

Figure 4. Model bias (model - observations) over depths using CTD data for short (2016–2018) experimental runs, highlighting temperature (a–c) and salinity (d–f) bias and CRMSE over depths (shading). HOTSSea v0.18 corresponds to the setup used in the final model run without bias correction (v1.01). HOTSSea v0.12 was omitted from the plots as results were nearly identical to v0.14. See Fig. 2 for a map of subdomains.

impacts of using coarse ERA5 atmospheric forcings. Among all experimental runs, v0.12 and v0.14 had the greatest nearsurface biases in temperature and salinity in the northern part of the model domain (i.e. DI subdomain; Fig. 4a, d) and in the central part of the domain (i.e. SGS domain; Fig. 4c, f). We investigated possible reasons and diagnosed that winds in ERA5 are generally weak and less variable throughout the domain relative to the higher-resolution HRDPS winds.
The HOTSSea v0.16 experiment helped isolate and evaluate the effect of using ORAS5 for boundary conditions at the Juan de Fuca Strait open-ocean boundary versus the relatively highly resolved CIOPS-W. In the northernmost subdomain, the Discovery Islands (DI), v0.16 showed similar biases to v0.14 with respect to salinity and temperature over depths (Fig. 4a, d). Looking to the SGS subdomain in the central part of the model domain, a relatively large warm bias appeared in v0.16 at depths greater than 10m (approximately +0.3 °C at 10m depth, increasing to +0.5 °C at 100+ m depths; Fig. 4c). In the JFS subdomain closest to the open-ocean boundary, the biases in temperature at depths >5 m were greater in v0.16 than in the other subdomains and greater relative to previous experimental runs, with biases of >1 °C near the surface and +0.4 to +0.75 °C at greater depths (Fig. 4b). The v0.16 experiment’s performance with respect to salinity in JFS was comparable to v0.14 (Fig. 4e). Our interpretation of these results is that in v0.16 the model inherited a warm bias in mid- and deep waters from the ORAS5 boundary conditions, evidenced by decreasing temperature bias with increasing distance from the open-ocean boundary. Contrary to expectations, using the lower-resolution ERA5 atmospheric forcings in the HOTSSea v0.14 experiment resulted in less of a warm temperature bias across all depths in the JFS subdomain at depths up to 100 m and had little discernible effect on salinity as compared with results of using the higher-resolution HRDPS outputs for atmospheric conditions in the v0.16 experiment. The better performance of v0.14 in JFS was unexpected because the model domain is relatively poorly resolved by ERA5 versus HRDPS. We suspect that wind- driven vertical mixing is biased low when using ERA5, and therefore ERA5 masks the effect of a near-surface warm bias introduced by ORAS5 by keeping the biased-warm water closer to the surface. Further investigation is warranted, and understanding the sources of biases is a priority.
In the final experimental run, HOTSSea v0.18, the RDRS v2.1 (∼ 10km horizontal) atmospheric conditions were used, representing an intermediate between the highest-resolution HRDPS forcings available only for 2014–2020 and the lower-resolution ERA5 forcings available for the entire hindcast period. The run performed similarly or better than the previous two runs with respect to temperature and salinity in most subdomains. The JFS subdomain was an exception, where the warm temperature bias was greater than in v0.14 (Fig. 4b). As above, the increased temperature bias was unexpected since v0.14 used the relatively weak and poorly resolved ERA5 atmospheric conditions. Thus, using both higher resolution and presumably more accurate atmospheric products (RDRS and HRDPS) resulted in greater temperature bias over depths up to ∼ 80 m in the JFS subdomain. This result is consistent with the diagnosis we reached earlier that the ERA5 atmospheric forcing fields used in the v0.14 experiment masked the effects of a warm temperature bias in the JFS subdomain up to depths of ∼ 10 m due to weaker winds and a shallower mixing as compared to the more powerful HRDPS and RDRS winds in the v0.14 and v0.18 experiments, respectively.
The preliminary experiments illuminated several important factors affecting model performance related to atmospheric and ocean boundary forcings. Coarsening the model bathymetry from 500 m to 1.5 km horizontally has especially impacted model performance in areas with narrow topography when compared with SalishSeaCast but otherwise has had only a minor effect on the evaluated aspects of model performance. Puget Sound (PS) remains essentially unevaluated (though see Tables S1 and S2) due to a lack of collated data for the area, but we suspect similar issues to the DI subdomain due to narrow topography; expanding our evaluation to PS remains a priority for future work. An important takeaway was an apparent systemic temperature bias present in ORAS5, the only reanalysis for hindcasting the ocean boundary conditions to 1980 and prior available at the time of writing. The ERA5 product similarly was deemed to be inadequate for the present purposes due to underpowered winds throughout the domain (RDRS v2.1 was used instead). The experiments prompted us to verify and quantify the bias introduced at the ocean boundary and to determine if simple bias correction could hold potential to reduce or eliminate inherited biases (results below).
4.2		Bias correction of open-ocean boundary conditions
The effect of the temperature bias corrections applied in HOTSSea v1.02 substantially improved model skill versus v1.01 (Figs. 5, 6; Tables S1, S2). In the JFS subdomain, for example, temperature bias calculated over all depths was reduced from +0.66 to +0.12 °C. The HOTSSea v1.01 run without bias correction had mean temperature biases calculated over all depths of +0.21 to +0.66 °C (excluding the PS subdomain with sparse observations) versus -0.14 to +0.2 °C in HOTSSea v1.02 (Table S1). Applying the temperature bias correction improved the model’s performance in the most distant subdomain from JFS, the DI subdomain – an unexpected result that emphasises the importance of accurate ocean boundary forcings.
4.3		Model–observation evaluation
4.3.1		Vertical profiles
Overall, the full hindcast (HOTSSea v1.02) performs best in the central and northern portion of the domain (Fig. 5). In the Strait of Georgia North (SGN) subdomain where the most CTD casts were available (N = 12 288), the temperature bias over all depths over the hindcast period was -0.08 °C (0–30 m = -0.39 °C; 30–150 m = -0.14 °C; >150 m = +0.13 °C) with an overall WSS of 0.97 and a correlation coefficient (R) of 0.94 (Tables S1, S2). Performance was similar in the SGS subdomain. The normalised target and Taylor diagrams (Figs. 6, S1–S3) indicate the model captures the seasonal ocean temperature variability well, though NCRMSE was >0.5 for JFS in the west and DI subdomains in the north. Model standard deviation is generally higher than observations in shallow depths and lower in deep water. Except for JFS and DI subdomains, the mean correlation coefficient with respect to temperature taken over all depths was >0.9 (Fig. S1a; Table S1). As depths increase, correlation generally decreases. At depths greater than 150 m, the model underestimates temperature variability, and the correlation coefficient was <0.8 for DI, JFS, and SGN. The southernmost PS subdomain was omitted from analysis due to relatively few collated observations, and preparing data for this area is a priority for future work. All metrics are tabulated by subdomain and depth strata in Tables S1 and S2 for reference.
HOTSSea v1.02 model skill was relatively poor with respect to salinity compared with temperature, especially in the DI and JFS subdomains (Table S2; Figs. S2–S3). However, performance was better in the Strait of Georgia, where mean bias taken over all depths in the SGS subdomain was +0.38 PSU (WSS = 0.98 and R = 0.96; Table S2). The narrow topography in the DI combined with the 1.5 km horizontal resolution likely leads to the observed error, as evidenced by experiments (Fig. 4). Another reason for salinity biases to increase in subdomains farther from the Fraser River could be that a climatology is used for estimating input from other rivers in the domain, whereas measurements are available for the Fraser River (Morrison et al., 2012). The open-ocean boundary conditions are forced using a climatology at the northern boundary in the DI subdomain, and this would affect model skill. However, the latter two explanations are considered less likely given the same boundary climatology and river forcings were used here as were used in SalishSeaCast, which performs well in the DI subdomain (Fig. 4).

Figure 5. Results of analysis using CTD data, aggregated by subdomain, showing mean temperature and salinity over depths (solid line: observations; dashed black line: HOTSSea v1.01 without bias correction; dashed red line: HOTSSea v1.02 with bias correction).

4.3.2		Sea surface temperature and salinity
Time series of SST and SSS taken at lightstations in the region are some of the longest in existence and present a valuable opportunity to evaluate model performance. SST and SSS were evaluated for the full hindcast using data from sampling at lightstations and buoys with temporal coverage of at least 10 years. The Taylor plots indicate the model reproduces the variability in SST well (Fig. 7a). Most stations have NCRMSE <0.5 with a correlation coefficient greater than 0.9. HOTSSea v1.02 performs well at stations within the Strait of Georgia and relatively poorly at Sheringham Point in the Juan de Fuca Strait, where the standard deviation of the model is approximately 20 % higher than observations. The normalised target diagram (Fig. 7c) shows that the model is typically biased cold with the mean bias typically <0.75 °C. The one lightstation with >1 °C bias was Chrome Island, which is located on an island close to shore where bathymetry is poorly resolved.
The evaluation of SSS indicates that HOTSSea v1.02 typically overestimates variability in SSS across the domain, as the normalised standard deviation of the model is greater than the observations at all stations except Chrome Island (Fig. 7b). Although the model’s performance with respect to SST is better than SSS, the model’s normalised standard deviation for SSS at all stations except Cape Mudge deviated by less than ±30 % of the observations. The correlation coefficient for SSS at most stations was relatively poor (0.4– 0.7) compared to SST (>0.85). The target plot indicates that the model is biased fresh at the surface by approximately -1 PSU at Cape Mudge in the northern part of the domain and -2 PSU at Entrance Island in the central Strait of Georgia, whereas at other stations the fresh bias was relatively small (< - 1 PSU; Fig. 7d).
Many applications of HOTSSea v1 in support of ecosystem modelling and research related to Pacific salmon are anticipated to require accuracy at weekly, monthly, or seasonal timescales rather than daily or hourly. To investigate whether the lower SSS performance was due to difficulty in capturing dynamics over hourly or daily timescales versus longer timescales, we applied a monthly moving average to both the model and observations, after which statistics were recalculated. Applying the moving average improved the evaluated statistics at all stations except Sheringham Point (Fig. 7, arrows). The improvement leads us to conclude that if research applications require bias <0.5 °C or <1.0 PSU, then averaging results to monthly timescales may ensure the accuracy is acceptable with respect to SST and SSS. As shown in Sect. 4.3.1, the model performance generally improves with depth.

Figure 6. Target plots of the model’s temperature bias and normalised centred root mean squared error (NCRMSE) using CTD data grouped by depth strata (panels a–d) and by subdomain. Results without bias correction to western open-ocean boundary conditions (HOTSSea v1.01) are shown in grey, and results with temperature bias correction (HOTSSea v1.02) are shown in red.

4.4		Decadal temperature trend evaluation
4.4.1		Nanoose station
The model represents seasonal changes in salinity and temperature over depths well. Both the modelled and observed biweekly climatologies at Nanoose station depict a characteristic intrusion of cold water and a temperature inversion that occurs in the area in the spring and late summer or fall (Fig. 8d, e), typically associated with upwelling events and neap tides (Johannessen et al., 2014; Masson, 2002; Riche et al., 2014). Strongly stratified and warm surface layers are shown in the shallower layers in the summer in both climatologies from the observations and from the model. The model–observation biases are generally largest at depths less than 3 m for both temperature and salinity. Although at depths >10 m there is only minor bias with respect to salinity (<0.5 PSU), salinity is biased fresh in depths <3 m, especially in the summer when it approaches -3 PSU (Fig. 8c).
The fresh bias observed here is not surprising; it was also apparent in the higher-resolution SalishSeaCast model (e.g. Fig. 4f; Soontiens and Allen, 2017) which formed the basis for HOTSSea; however, it remains an important area for improvement in the future and may affect model skill with respect to circulation patterns. The modelled temperatures are biased warm in the summer and cool in the winter by a maximum of approximately 1 °C at depths less than 3 m. In contrast, the modelled temperatures are consistently biased slightly high (∼ 0.5 °C) at depths over 100 m.
After truncating the Nanoose station time series to match the model hindcast period, the monotonic temperature trend across depths 4–400 m was calculated using the MK test and the TS slope estimator. Significant autocorrelation was detected in the detrended biweekly time series, so the 3PW algorithm was used to adjust for this (see Sect. 3.4). The trend from observations at Nanoose station was evaluated to be 0.031 °C per decade but was not found to be statistically significant (LCL: -0.018; UCL: 0.08; p ≤ 0.12). HOTSSea

Figure 7. Taylor (a, b) and target (c, d) plots evaluating HOTSSea v1.02 sea surface temperature (SST; a, c) and salinity (SSS; b, d). SST was evaluated using data from buoys (numbered) and lightstations (named), and SSS was evaluated using lightstation data only. The standard deviation (solid grey contours in Taylor diagrams) and centred root mean square error (dashed grey contours in Taylor diagrams) have been normalised to enable comparison on the same plot. Perfect agreement between model and observations on the Taylor plot would correspond to normalised standard deviation = 1, correlation = 1, and NCRMSE = 0. Perfect agreement between model and observations on the target plot would be at the centre. Arrows depict the change after applying a 1-month moving average on both the model and observation time series.

v1.01 without ocean boundary temperature bias correction detected no trend (0.00 °C per decade; LCL: -0.06; UCL: 0.06; p ≤ 0.73). HOTSSea v1.02 with ocean temperature bias correction predicted an insignificant trend of 0.012 °C per decade (LCL: -0.052; UCL: 0.075; p ≤ 0.3). Although trends over the water column for the 1980–2018 hindcast period were not significant at the 95 % confidence level, results suggested that HOTSSea v1.02 with ocean boundary temperature bias correction performs better than HOTSSea v1.01 with respect to long-term trends.
The observed inter-annual and intra-annual variability over the hindcast period is well captured by the model (Fig. 9). The largest deviations between modelled (HOTSSea v1.02) and observed temperature anomalies over the whole water column were <0.5 °C. Relatively large deviations from observations in the 1980–1983 period suggest that the 1-year spin-up may be too short, and this remains a priority area to investigate in the future. Warm anomalies as a result of a documented oceanic heat wave circa 2015–2016 (Gruber et al., 2021; Khangaonkar et al., 2021b) are generally well represented, suggesting the model is capable of reproducing extreme events. The evaluation of the trends at Nanoose station also indicate the model does not incur a detectable drift in temperature bias at this location, despite being run with no data assimilation.
The secular trend previously calculated for the 1970 to 2005 period of 0.24 °C per decade (95% CI ± 0.1 °C; Masson and Cummins, 2007) was over 6-fold greater than calculated here for the 1980 to 2018 hindcast period. We verified that the difference was not attributable to the methods (TS vs. LR for slope estimation) by re-estimating the trend from the observations using LR for the same 1980– 2018 period, finding a similar result using LR (0.038 °C per decade). The weaker secular trend calculated for the 1980–

Figure 8. Biweekly climatologies and model–observation bias over depths over the hindcast period using measurements taken at Nanoose station between 1980 and 2018 (N = 5692). The four panels to the left (a, b, d, e) show salinity extracted from observations (a, d) and the model (b, e). Temperatures are similarly represented from observations and the model in the bottom plots. The two plots on the right show the model bias (i.e. model - observations) over depth. Data were binned bimonthly.

Figure 9. Temperature anomalies (seasonal) from observations (black, solid) at Nanoose station in the central Strait of Georgia versus those derived from HOTSSea v1.02 model outputs (red, dashed), depth integrated over 4.5–400 m. The grey area represents the model bias. Observations from the 1970–1980 period at Nanoose station are included to illustrate the cooler period, a multiyear swing in temperature anomalies (1977, 1978), followed by a regime shift occurring circa 1977 (Beamish et al., 1999; Hare and Mantua, 2000).

2018 period is in part due to the omission of the 1970s, a period with cold water temperature anomalies observed at Nanoose station (Figs. 9, S4). The period of the late 1970s corresponds approximately to a polarity shift in the Pacific Decadal Oscillation (PDO; Mantua et al., 1997; Mantua and Hare, 2002), a trough in the North Pacific Gyre Oscillation (NPGO; Di Lorenzo et al., 2008), which oscillates on decadal or multi-decadal timescales. The period of the late 1970s also corresponds to a trough in a multi-regime index circa 1977 that correlates with a reduction in northern Pacific fisheries catches (Beamish et al., 1999). Overall, the analysis reveals a key limitation of HOTSSea v1 – omission of the 1970s – which is thus a priority area for subsequent iterations. Our analysis also serves as an update to the MC07 analysis and confirms ocean temperatures have not subsequently reverted.
Analysis of trends at each depth at Nanoose station indicates that the modelled trends are in generally good agreement with observed trends (Fig. 10). Both the modelled and observed seasonal trends in the upper water column in the summer, fall, and annually were statistically significant, whereas the trends in the winter and spring were not. A deviation between the modelled and observed trends is apparent, especially in the summer, between approximately 20 and 100 m. The reason for the deviation is yet to be determined, but a similar deviation is also apparent in the analysis using CTDs in the Gulf Islands (GI) subdomain (Fig. 5). The seasonal trend patterns over depth at Nanoose station indicate that trends vary seasonally and spatially (over depths) at this location.
4.4.2		Strait of Georgia
Modelled trends from HOTSSea v1.02 were evaluated in each grid cell in the central and northern part of the model domain (Strait of Georgia and surrounding waters) over the 1980–2018 hindcast period (Fig. 11). The 0–30 m depth stratum was selected given that the evaluation of the modelled versus observed trends showed generally good agreement in approximately the top 30m of the water column. This stratum also has special relevance for phytoplankton bloom dynamics (Allen and Wolfe, 2013). The model indicates the top 30 m is generally warming in the Strait of Georgia, with some areas experiencing statistically significant warming over the entire water column. Despite finding that the mean trend computed across all depths at Nanoose station in the winter was not statistically significant (Figs. 10a, 11f), a statistically significant warming trend in the winter was detected in the 0– 30 m depth stratum in the northeastern Strait of Georgia and in Jervis Inlet (Fig. 11a). The rapid warming of Jervis Inlet is a finding consistent with trends reported for other fjords in the region (Jackson et al., 2021). The analysis of temperature trends in the spring revealed no significant warming trend (Fig. 11b), similar to the analysis of the Nanoose data (Fig. 10b). In the summer, only the waters in the central and northern Strait of Georgia generally showed a statistically significant warming trend, whereas the southern Gulf islands and surrounding waters showed no statistically significant trend (Fig. 11c). The fall season may also have experienced the most spatially consistent secular ocean temperature warming since 1980 over the 0–30 m depth range - significant trends in the fall were generally widespread apart from Howe Sound (Fig. 11d). Annualised trends were relatively low but generally statistically significant throughout the Strait of Georgia (Fig. 11e). Note that the magnitudes of all modelled trends throughout the domain are likely biased low, as the analysis of the modelled secular trend at Nanoose station was weak relative to the observed secular trend (see above).
5		Discussion and conclusion
To our knowledge, HOTSSea v1 represents the first 3D physical oceanographic model for the Salish Sea to be extended back to 1980, as well as the first model to use recently available forcing products to do so (e.g. RDRS v2.1, ORAS5). A highlight of the results of our evaluation of HOTSSea v1 was the skill with which the model represents observed decadal- scale ocean temperature trends at Nanoose station – without incurring any detectable drift despite being a free run with no data assimilation (Figs. 9, 10). After confirming the model’s skill at recreating the trend at Nanoose station, we used the model to further our understanding of trends in areas within the central and northern part of the domain with sparse observations, offering a first glimpse at the spatial– temporal heterogeneity of ocean temperature trends in the area (Fig. 11). The 0–30 m depth range was of particular interest because changes occurring at these depths affect dynamics of spring and fall phytoplankton blooms (Masson and Peña, 2009). Most importantly, HOTSSea v1 outputs indicate that the warming trend apparent at Nanoose station is not an isolated phenomenon within the Salish Sea. Gradual ocean temperature changes such as those quantified by our analyses have the potential to affect the growth, body mass, and marine distributions of fish (Pauly, 2021; Pauly and Cheung, 2018) and support research related to the velocity of climate change (Burrows et al., 2014; Loarie et al., 2009). Anomalous events such as marine heat waves are expected to occur more frequently due to climate change, and these events can lead to key biophysical thresholds being surpassed, impacting marine organisms (Gruber et al., 2021). Given that HOTSSea v1 generally replicates observed ocean temperature anomalies and extremes, it therefore shows promise for investigating climate-related pathways of effects from physical ocean properties to marine ecosystems.
Overall, our evaluation quantifies the biases of the model at several spatial–temporal scales, with special attention to those scales and variables anticipated to be relevant for marine ecosystem model development, marine ecosystem research, and ecosystem-based fisheries management. Although variability in physical conditions at fine spatial– temporal scales (i.e. minutes to days, metres to kilometres) is an important determinant of the structure and dynamics of marine ecosystems (e.g. Jones et al., 2014), many research applications related to marine ecosystems focus on processes and mechanisms that exhibit variability on weekly or greater timescales and at spatial scales of tens of kilometres (Fulton et al., 2019). Locally, juvenile salmon enter the marine environment from natal rivers over a protracted period of weeks to months (Groot and Margolis, 1991; Healey, 1980); migrate at rates of 1.5–20 km d-1 (Melnychuk et al., 2010,

Figure 10. Seasonal and annual trends (1980–2018) over depths observed at Nanoose station compared with HOTSSea v1.02 outputs.
Shading represents upper and lower 95 % confidence limits.

Figure 11. Mean ocean temperature trends in the Strait of Georgia computed as seasonal mean trends (columns) for two depth strata (rows), 0–30 m depth stratum and over all depths (all z), from HOTSSea v1.02 model outputs. A mask (grey) has been applied to grid cells that are shallower than the depth stratum. Grey hatching has been applied to grid cells where the trend was not statistically significant. The star symbol denotes the approximate location of Nanoose station.
2013; Trudel et al., 2009; Tucker et al., 2009; Welch et al., 2011); are widely distributed throughout the pelagic zone (Beamish et al., 2011; Riddell et al., 2018; Trudel and Tucker, 2013); and exhibit weekly, monthly, and seasonal variability in growth and mortality. We therefore expect that hindcasting physical ocean water properties at the spatial–temporal resolution used here will be of value for many studies related to Pacific salmon.
A contemporary challenge to achieving a long hindcast for the Salish Sea is related to the quality of the available products to use for boundary and atmospheric forcing. The experimental approach we used helped us isolate the effect of forcings and our model’s horizontal resolution from other factors affecting model performance, revealing substantial biases inherited from both atmospheric and oceanic boundary forcings. What we learned using this approach will guide future efforts to improve and extend HOTSSea v1 and possibly other models for the domain. Using the RDRS v2.1 forcings in the final hindcast led to substantially better model performance versus using ERA5 (Sect. 4.1) – presently the only alternative to our knowledge with coverage of the full hindcast period. Wind speed biases in ERA5 have been noted in other coastal and mountainous areas (Potisomporn et al., 2023). The results of our experiments similarly indicated that locally weak winds in ERA5 resulted in reduced mixing and increased near-surface biases, leading to especially poor model performance in areas with complex and narrow topography (e.g. the Discovery Islands). Winds are an important factor determining total productivity of the Strait of Georgia (Collins et al., 2009; Johannessen et al., 2020), and thus it was fortunate that the RDRS v2.1 atmospheric reanalysis was recently made available. The experiments also revealed that the global ocean reanalysis ORAS5 outputs exhibit temperature and salinity biases at the mouth of the Juan de Fuca Strait which were inherited to some extent by our model. We suspect the main issue is that the Salish Sea is poorly resolved in ORAS5, leading to poor representation of estuarine flow and physical dynamics at the mouth of the Juan de Fuca Strait. The temperature bias correction we applied here to ORAS5 at the Juan de Fuca Strait boundary improved model performance substantially, even in areas of the domain far from the ocean boundary (Figs. 5, 6). The success of our crude bias correction serves to highlight the potential benefit of relatively advanced techniques such as statistical downscaling using machine learning or the application of data assimilation techniques to produce a reanalysis (Adachi and Tomita, 2020). Another avenue worth exploring would be to use outputs from other regional ocean hindcasts as boundary forcings, should they become available (e.g. Lin et al., 2022; Peña et al., 2016). Applying bias correction to salinity at the boundary could also be of benefit and may help with improving the model’s near-surface fresh bias.
Next steps include several improvements, extensions, and applications. Our preliminary evaluation of the model prioritised the Juan de Fuca Strait and the central part of the Salish Sea (i.e. the Strait of Georgia and surrounding waters), whereas Puget Sound remains essentially unevaluated (though see Tables S1 and S2). Collated Canadian observations had relatively sparse coverage of Puget Sound (Fig. 2, Table 2), and adding more observations from American institutions therefore is a near-term priority. The reported biases in the final model are also a priority to diagnose and correct, especially near-surface fresh biases in summer months (e.g. Fig. 8c). Similar near-surface biases were reported in Salish- SeaCast (Soontiens and Allen, 2017). The salinity bias is of particular concern as it may hinder the model’s capability of representing circulation patterns. Evaluation of tides and the estuarine circulation is a particularly important task, as these are primary factors affecting deepwater renewal in the Sal- ish Sea (Ebbesmeyer et al., 1989; MacCready et al., 2021). Once circulation is evaluated, one potentially promising application would be to use velocity fields from HOTSSea v1 for Lagrangian particle simulations (Hernández-Carrasco et al., 2020; Snauffer et al., 2014). Changes to regional-scale atmospheric patterns may have affected sea surface height and circulation patterns, which in turn may affect important ecosystem processes such as larval dispersal. A logical next step is to use data assimilation to combine information from dynamical simulation with observations to produce a reanalyses, with estimates for the oceanic fields that are maximally consistent with observations (Zaron, 2011). Data assimilation techniques are developing rapidly and show great promise for helping minimise model error. Finally, upgrading from NEMO v3.6 to NEMO v4.x is a priority, as new vertical mixing options and a wetting-and-drying scheme for intertidal areas hold potential for improving model skill in this domain (Madec et al., 2023).
Ultimately, we hope to extend HOTSSea further back in time, motivated by a need to capture the period prior to an oceanic regime shift that occurred circa 1977 that can be seen in the Nanoose station time series (Figs. 9, S4). The regime shift is of particular interest for fisheries-related research in the area (Beamish et al., 1999; Di Lorenzo et al., 2008; Mantua et al., 1997; Mantua and Hare, 2002; MC07). Despite various climate regime indices subsequently reverting in polarity, ocean temperatures may not have reverted. Coupling the physical model with a biogeochemical module is another potentially valuable avenue to explore, as integrating biogeochemistry into the hindcast would provide retrospective details of oxygen, nutrients, carbonate chemistry, aragonite saturation, and pH – all of which have changed since pre-industrial times in the Salish Sea (Jarníková et al., 2022). Secular ocean warming and increased seasonal variability have also had a quantifiable effect on oxygenation of deep fjord waters in the region (Jackson et al., 2021a), which may have affected the productivity or community composition of lower trophic levels (Johannessen et al., 2020). However, computational costs associated with traditional biogeochemical models are currently a challenge that would need to be overcome. The HOTSSea v1 reconstruction of temperature fields over depths should also prove useful for investigating the relative impact of pathogens and predation on fish. Warming ocean temperatures lead to greater vulnerability to pathogens and disease and thus to greater vulnerability to predation (Miller et al., 2014; Teffer et al., 2018), but concurrent increases in predator abundances make it difficult to determine whether predation is a primary or secondary mortality factor (Walters and Christensen, 2019). Available observational time series are mainly limited to surface waters, and HOTSSea v1 can help address this gap by providing ocean temperature fields at specific depths and areas occupied by fish.
In conclusion, HOTSSea v1 shows promise for improving our understanding of long-term physical changes occurring in the Salish Sea with relevance for fish, fisheries, and marine ecosystems. The model’s performance with respect to salinity, temperature, decadal-scale trends, and temperature anomalies indicates HOTSSea v1 can support ecosystem research focused on dynamics that unfold over months, seasons, years, or decades. The model addresses significant gaps in historical observations and can help drive biogeochemical and ecosystem models aimed at revealing dominant drivers of ecological productivity (Hermann et al., 2023; Macias et al., 2014; Piroddi et al., 2021).
Code and data availability. HOTSSea is based on the NEMO source code version 3.6 (https://doi.org/10.5281/zenodo.1464816, Madec et al., 2017; subversion trunk revision rev10584), released under the open-source CeCILL license (https://cecill.info, last access: 7 March 2024). The HOTSSea v1 configuration files and source code used for analysis and visuals in the present article have been archived at https://doi.org/10.5281/zenodo.13887813 (Oldford, 2024a). Forcings and observations prepared for the model are provided separately due to space limits, with HRDPS at https: //doi.org/10.5281/zenodo.12193924 (Oldford and Dunphy, 2024a); RDRS at https://doi.org/10.5281/zenodo.12206291 (ECCC, 2024); and ERA5, CIOPS, runoff, and observational data at https://doi.org/ 10.5281/zenodo.12312769 (Oldford and Dunphy, 2024b). Raw outputs are currently stored on a server without necessary bandwidth to make them publicly available given the file sizes; however, they can be made available by contacting the authors. Depth-averaged, monthly mean outputs are available online at https://doi.org/10. 5281/zenodo.13942109 (Oldford, 2024b). The processed monthly HOTSSea outputs may be accessed by installing the pacea opensource R package, which aims to reduce the burden of wrangling such data for fisheries stock assessments and other applications (https://doi.org/10.5281/zenodo.13840804, Edwards et al., 2024).
Supplement. The supplement related to this article is available online at: https://doi.org/10.5194/gmd-18-211-2025-supplement.
Author contributions. All authors contributed to the review and editing of the paper. GO wrote the original draft and conducted the data preparation, analysis, formal evaluation, and data visualisations. MD did the server configuration, compilation, and running of the model with contributions from GO. TJ provided expertise and carried out analyses related to trend analysis and climatologies. VC and MD contributed to the administration and supervision of the model development.
Competing interests. The contact author has declared that none of the authors has any competing interests.
Disclaimer. Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
Acknowledgements. The authors thank Susan Allen, Andrew Edwards, and Carl Walters for their review and suggestions on drafts of the manuscript and extend thanks to Elise Olson for early discussions about analysis and code. The authors also acknowledge the work of Susan Allen and colleagues, who developed the SalishSea- Cast model for the same domain which served as the starting point for the HOTSSea v1 model.
Financial support. This research has been supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) Mitacs programme (grant no. IT12313), Fisheries and Oceans Canada (Government of Canada), the Pacific Salmon Foundation, the University of British Columbia (UBC), the UBC Ocean Leaders programme, and the Salish Sea Marine Survival Project.
Review statement. This paper was edited by Vassilios Vervatis and reviewed by two anonymous referees.
References
Adachi, S. A. and Tomita, H.: Methodology of the Constraint Condition in Dynamical Downscaling for Regional Climate Evaluation: A Review, J. Geophys. Res.-Atmos., 125, e2019JD032166, https://doi.org/10.1029/2019JD032166, 2020.
Allen, S. E. and Wolfe, M. A.: Hindcast of the timing of the spring phytoplankton bloom in the strait of Georgia, 1968–2010, Prog. Oceanogr., 115, 6–13, https://doi.org/10.1016/j.pocean.2013.05.026, 2013.
Amos, C. L., Martino, S., Sutherland, T. F., and Al Rashidi, T.: Sea Surface Temperature Trends in the Coastal Zone of British Columbia, Canada, J. Coast. Res., 300, 434–446, https://doi.org/10.2112/JCOASTRES-D-14-00114.1, 2015.
Araujo, H. A., Holt, C., Curtis, J. M. R., Perry, R. I., Irvine, J. R., and Michielsens, C. G. J.: Building an ecosystem model using mismatched and fragmented data: a probabilistic network of early marine survival for coho salmon Oncorhynchus kisutch in the Strait of Georgia, Prog. Oceanogr., 115, 41–52, https://doi.org/10.1016/j.pocean.2013.05.022, 2013.
Arakawa, A. and Lamb, V. R.: Computational design of the basic dynamical processes of the UCLA general circulation model, Gen. Circ. Models Atmosphere, 17, 173–265, 1977.
Beamish, R. J.: Climate and Exceptional Fish Production off the West Coast of North America, Can. J. Fish. Aquat. Sci., 50, 2270–2291, https://doi.org/10.1139/f93-252, 1993.
Beamish, R. J. (Ed.): Climate Change and Northern Fish Populations, National Research Council of Canada, ISBN: 0-66015780-2, Ottawa, ON, Canada, 739 pp., 1995.
Beamish, R. J., Noakes, D. J., McFarlane, G. A., Klyashtorin, L., Ivanov, V. V., and Kurashov, V.: The regime concept and natural trends in the production of Pacific salmon, Can. J. Fish. Aquat. Sci., 56, 516–526, https://doi.org/10.1139/f98-200, 1999.
Beamish, R. J., Sweeting, R. M., Lange, K. L., Noakes, D. J., Preikshot, D. B., and Neville, C. M.: Early Marine Survival of Coho Salmon in the Strait of Georgia Declines to Very Low Levels, Mar. Coast. Fish., 2, 424–439, https://doi.org/10.1577/C09- 040.1, 2010.
Beamish, R. J., Lange, K. L., Neville, C.-E. M., Sweeting, R. M., and Beacham, T. D.: Structural patterns in the distribution of ocean- and stream-type juvenile chinook salmon populations in the Strait of Georgia in 2010 during the critical early marine period, North Pacific Anadromous Fish Commission, https://www.npafc.org/wp-content/uploads/ Public-Documents/2011/1354Canada.pdf (last access: 14 January 2025), 2011.
Bluestein, L.: A linear filtering approach to the computation of discrete Fourier transform, IEEE Trans. Audio Electroacoustics, 18, 451–455, https://doi.org/10.1109/TAU.1970.1162132, 1970.
Brockwell, P. J. and Davis, R. A.: Introduction to Time Series and Forecasting, Springer International Publishing, Cham, https://doi.org/10.1007/978-3-319-29854-2, 2016.
Burrows, M. T., Schoeman, D. S., Richardson, A. J., Molinos, J. G., Hoffmann, A., Buckley, L. B., Moore, P. J., Brown, C. J., Bruno, J. F., Duarte, C. M., Halpern, B. S., Hoegh-Guldberg, O., Kappel, C. V., Kiessling, W., O’Connor, M. I., Pandolfi, J. M., Parmesan, C., Sydeman, W. J., Ferrier, S., Williams, K. J., and Poloczanska, E. S.: Geographical limits to species-range shifts are suggested by climate velocity, Nature, 507, 492–495, https://doi.org/10.1038/nature12976, 2014.
Collaud Coen, M., Andrews, E., Bigi, A., Martucci, G., Roma- nens, G., Vogt, F. P. A., and Vuilleumier, L.: Effects of the prewhitening method, the time granularity, and the time segmentation on the Mann–Kendall trend detection and the associated Sen’s slope, Atmos. Meas. Tech., 13, 6945–6964, https://doi.org/10.5194/amt-13-6945-2020, 2020.
Collins, A. K., Allen, S. E., Pawlowicz, R., Collins, A. K., Allen, S. E., and Pawlowicz, R.: The role of wind in determining the timing of the spring bloom in the Strait of Georgia, Can. J. Fish. Aquat. Sci., 66, 1597–1616, https://doi.org/10.1139/F09- 071, 2009.
Cooley, J. W. and Tukey, J. W.: An Algorithm for the Machine Calculation of Complex Fourier Series, Math. Comput., 19, 297– 301, 1965.
Copernicus Climate Change Service: ORAS5 global ocean reanalysis monthly data from 1958 to present, Copernicus Climate Change Service (C3S) Climate Data Store (CDS) [data set], https://doi.org/10.24381/CDS.67E8EEB7, 2021.
Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M. A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A. C. M., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A. J., Haimberger, L., Healy, S. B., Hersbach, H., Hólm, E. V., Isaksen, L., Kållberg, P., Köhler, M., Matricardi, M., Mcnally, A. P., Monge-Sanz, B. M., Morcrette, J. J., Park, B. K., Peubey, C., de Rosnay, P., Tavolato, C., Thépaut, J. N., and Vitart, F.: The ERA-Interim reanalysis: Configuration and performance of the data assimilation system, Q. J. Roy. Meteor. Soc., 137, 553–597, https://doi.org/10.1002/qj.828, 2011.
de Mutsert, K., Coll, M., Steenbeek, J., Ainsworth, C., Buszowski, J., Chagaris, D., Christensen, V., Heymans, S. J. J., Lewis, K. A., Libralato, S., Oldford, G., Piroddi, C., Romagnoni, G., Serpetti, N., Spence, M. A., and Walters, C.: Advances in spatial-temporal coastal and marine ecosystem modeling using Ecospace, in: Reference Module in Earth Systems and Environmental Sciences, Elsevier, https://doi.org/10.1016/B978-0-323-90798-9.00035-4, 2023.
Digital Research Alliance of Canada: Advanced Research Computing, https://alliancecan.ca/en/services/ advanced-research-computing, last access: 10 June 2022.
Di Lorenzo, E., Schneider, N., Cobb, K. M., Franks, P. J. S., Chhak, K., Miller, A. J., McWilliams, J. C., Bograd, S. J., Arango, H., Curchitser, E., Powell, T. M., and Rivière, P.: North Pacific Gyre Oscillation links ocean climate and ecosystem change, Geophys. Res. Lett., 35, L08607, https://doi.org/10.1029/2007GL032838, 2008.
Dosser, H. V., Waterman, S., Jackson, J. M., Hannah, C., and Hunt, B. P. V.: Tidal mixing maintains regional differences in water properties and nutrient ratios in British Columbia coastal waters, ESS Open Archive [preprint], https://doi.org/10.1002/essoar.10505244.1, 2020.
Dosser, H. V., Waterman, S., Jackson, J. M., Hannah, C. G., Evans, W., and Hunt, B. P. V.: Stark Physical and Biogeochemical Differences and Implications for Ecosystem Stressors in the Northeast Pacific Coastal Ocean, J. Geophys. Res.-Oceans, 126, e2020JC017033, https://doi.org/10.1029/2020JC017033, 2021.
Ebbesmeyer, C. C., Coomes, C. A., Cannon, G. A., and Bretschnei- der, D. E.: Linkage of Ocean and Fjord Dynamics at Decadal Period, in: Aspects of Climate Variability in the Pacific and Western Americas, edited by: Peterson, D. H., American Geophysical Union (AGU), 399–417, https://doi.org/10.1029/gm055p0399, 1989.
Edwards, A., Tai, Travis, Watson, J., Peña, M. A., Hilborn, A., Hannah, C. G., Rooper, C. N., Flynn, K. L., and Oldford, G.: pacea: An R package of Pacific ecosystem information to help facilitate an ecosystem approach to fisheries management, v1.0.0, Zenodo [data set], https://doi.org/10.5281/zenodo.13840804, 2024.
Environment and Climate Change Canada (ECCC): High Resolution Deterministic Prediction System (HRDPS) – Continental [data set] (7.0.0), ECCC Identifier: 5b401fa0-6c29- 57f0-b3d5-749f301d829d, https://eccc-msc.github.io/open-data/ msc-data/nwp_hrdps/readme_hrdps-datamart_en/ (last access: 10 January 2023), 2020.
Environment and Climate Change Canada (ECCC): HOTSSea Forcings – RDRS v2.1, Zenodo [data set], https://doi.org/10.5281/zenodo.12206291, 2024.
Esenkulova, S., Suchy, K. D., Pawlowicz, R., Costa, M., and Pearsall, I. A.: Harmful Algae and Oceanographic Conditions in the Strait of Georgia, Canada Based on Citizen Science Monitoring, Front. Mar. Sci., 8, 1193, https://doi.org/10.3389/FMARS.2021.725092, 2021.
Fatland, R., MacCready, P., and Oscar, N.: Chapter 14 – LiveO- cean, in: Cloud Computing in Ocean and Atmospheric Sciences, edited by: Vance, T. C., Merati, N., Yang, C., and Yuan, M., Academic Press, 277–296, https://doi.org/10.1016/B978-0-12- 803192-6.00014-1, 2016.
Feely, R., Doney, S., and Cooley, S.: Ocean Acidification: Present Conditions and Future Changes in a High-CO2 World, Oceanography, 22, 36–47, https://doi.org/10.5670/oceanog.2009.95, 2009.
Fisheries and Oceans Canada (DFO): British Columbia lightstation daily sea-surface temperature and salinity measurements, 1914 to present, DFO [data set], https://catalogue.cioospacific.ca/dataset/ ca-cioos_654a4ece-7271-4f5a-ba60-b50a62dbd051?local=en (last access: 10 January 2025), 2022a.
Fisheries and Oceans Canada (DFO): Vertical profiles of seawater properties measured by Conductivity-Temperature- Depth loggers in British Columbia, Canada, 1965 to present, https://catalogue.cioospacific.ca/dataset/ca-cioos_ 89aa45fc-7b42-426e-9293-e5703534bc4f?local=en (last access: 14 December 2023), 2022b.
Fisheries and Oceans Canada (DFO): Marine Environmental Data Section Archive, Wave Data Online, DFO [data set], https://meds- sdmm.dfo- mpo.gc.ca/isdm- gdsi/waves- vagues/ data-donnees/index-eng.asp (last access: 20 January 2023), 2024.
Foreman, M. G. G., Crawford, W. R., Cherniawsky, J. Y., Henry, R. F., and Tarbotton, M. R.: A high-resolution assimilating tidal model for the northeast Pacific Ocean, J. Geophys. Res.-Oceans, 105,	28629–28651,	https://doi.org/10.1029/1999JC000122, 2000.
Fulton, E. A., Blanchard, J. L., Melbourne-Thomas, J., Plagányi, É. E., and Tulloch, J. D. V.: Where the Ecological Gaps Remain, a Modelers’ Perspective, Front. Ecol. Evol., 7, 424, https://doi.org/10.3389/fevo.2019.00424, 2019.
Georgia Strait Alliance: About the Strait, https://georgiastrait.org/ issues/about- the- strait- 2/, last access: 17 January 2020.
Gasset, N., Fortin, V., Dimitrijevic, M., Carrera, M., Bilodeau, B., Muncaster, R., Gaborit, É., Roy, G., Pentcheva, N., Bulat, M., Wang, X., Pavlovic, R., Lespinas, F., Khedhaouiria, D., and Mai, J.: A 10km North American precipitation and land-surface reanalysis based on the GEM atmospheric model, Hydrol. Earth Syst. Sci., 25, 4917–4945, https://doi.org/10.5194/hess-25-4917- 2021, 2021.
Gocic, M. and Trajkovic, S.: Analysis of changes in meteorological variables using Mann-Kendall and Sen’s slope estimator statistical tests in Serbia, Glob. Planet. Change, 100, 172–182, https://doi.org/10.1016/j.gloplacha.2012.10.014, 2013.
Goldfeld, S. M. and Quandt, R. E.: Some Tests for Ho- moscedasticity,	J.	Am. Stat. Assoc., 60,	539–547, https://doi.org/10.1080/01621459.1965.10480811, 1965.
Gower, J., King, S., Statham, S., Fox, R., and Young, E.: The malaspina dragon: A newly-discovered pattern of the early spring bloom in the strait of georgia, british columbia, canada, Prog. Oceanogr., 115, 181–188, https://doi.org/10.1016/j.pocean.2013.05.024, 2013.
Groot, C. and Margolis, L.: Pacific salmon life histories, UBC press, Vancouver, BC, xv + 564 pp., ISBN 0-7748-0359-2, 1991.
Gruber, N., Boyd, P. W., Frölicher, T. L., and Vogt, M.: Biogeochemical extremes and compound events in the ocean, Nature, 600, 395–407, https://doi.org/10.1038/s41586-021-03981- 7,	2021.
Guan, L.: Spatial and temporal dynamics of larval fish assemblages in the Strait of Georgia, PhD Dissertation, University of Victoria, Victoria, BC, Canada, xvi + 158 pp., https://dspace.library.uvic. ca/items/72db1611-a9db-4d0e-920c-3b3ce854fd8f (last access: 10 January 2025), 2015.
Hakai Institute: Hakai Water Properties Vertical Profile Data Measured by Oceanographic Profilers, Provisional, Hakai Institute [data set], https://catalogue.hakai.org/dataset/ca-cioos_ 6143028b-028d-46c7-a67d-f3a513435e63?local=en (last access: 10 September 2024), 2024.
Hare, S. R. and Mantua, N. J.: Empirical evidence for North Pacific regime shifts in 1977 and 1989, Prog. Oceanogr., 47, 103–145, https://doi.org/10.1016/S0079-6611(00)00033-1, 2000.
Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N. J., Kern, R., Picus, M., Hoyer, S., van Kerkwijk, M. H., Brett, M., Haldane, A., del Río, J. F., Wiebe, M., Peterson, P., Gérard- Marchant, P., Sheppard, K., Reddy, T., Weckesser, W., Abbasi, H., Gohlke, C., and Oliphant, T. E.: Array programming with NumPy, Nature, 585, 357–362, https://doi.org/10.1038/s41586- 020-2649-2, 2020.
Harrison, P. J., Fulton, J. D., Taylor, F. J. R., Parsons, T. R., Futon, J. E., and Taylor, F. J. R.: Review of the Biological Oceanography of the Strait of Georgia: Pelagic Environment, Can. J. Fish. Aquat. Sci., 40, 1064–1094, https://doi.org/10.1139/f83- 129, 1983.
Healey, M. C.: The ecology of juvenile salmon in Georgia Strait, British Columbia, in: Salmonid ecosystems of the North Pacific, 203–229, ISBN 0870713353, 1980.
Helsel, D. R. and Hirsch, R. M.: Statistical methods in water resources, Elsevier, Netherlands, 522 pp., ISBN 0-444-81463-9, 1992.
Hermann, A. J., Cheng, W., Stabeno, P. J., Pilcher, D. J., Kearney, K. A., and Holsman, K. K.: Applications of Biophysical Modeling to Pacific High-Latitude Ecosystems, Oceanography, 36, 101– 108, 2023.
Hernández-Carrasco, I., Alou-Font, E., Dumont, P.-A., Cabornero, A., Allen, J., and Orfila, A.: Lagrangian flow effects on phytoplankton abundance and composition along filament-like structures, Prog. Oceanogr., 189, 102469, https://doi.org/10.1016/j.pocean.2020.102469, 2020.
Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R. J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., de Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., and Thépaut, J. N.: The ERA5 global reanalysis, Q. J. Roy. Meteor. Soc., 146, 1999– 2049, https://doi.org/10.1002/qj.3803, 2020.
Ianson, D., Allen, S. E., Moore-Maley, B. L., Johannessen, S. C., and Macdonald, and Robie W.: Vulnerability of a semienclosed estuarine sea to ocean acidification in contrast with hypoxia, Geophys. Res. Lett., 43, 5793–5801, https://doi.org/10.1002/2016GL068996, 2016.
Institute for Ocean Sciences (IOS) Water Properties Group, Water Properties Database, https://www.waterproperties.ca/, last access: 10 January 2025.
Islam, S. U., Hay, R. W., Déry, S. J., and Booth, B. P.: Modelling the impacts of climate change on riverine thermal regimes in western Canada’s largest Pacific watershed, Sci. Rep., 9, 11398, https://doi.org/10.1038/s41598-019-47804-2, 2019.
Jackson, J. M., Bianucci, L., Hannah, C. G., Carmack, E. C., and Barrette, J.: Deep Waters in British Columbia Mainland Fjords Show Rapid Warming and Deoxygenation From 1951 to 2020, Geophys. Res. Lett., 48, e2020GL091094, https://doi.org/10.1029/2020GL091094, 2021.
Jarníková, T., Ianson, D., Allen, S. E., Shao, A. E., and Olson, E. M.: Anthropogenic Carbon Increase has Caused Critical Shifts in Aragonite Saturation Across a Sensitive Coastal System, Glob. Biogeochem. Cycles, 36, e2021GB007024, https://doi.org/10.1029/2021GB007024, 2022.
Johannessen, S. C., Masson, D., and Macdonald, R. W.: Oxygen in the deep Strait of Georgia, 1951–2009: The roles of mixing, deep-water renewal, and remineralization of organic carbon, Limnol. Oceanogr., 59, 211–222, https://doi.org/10.4319/lo.2014.59.1.0211, 2014.
Johannessen, S. C., Macdonald, R. W., and Strivens, J. E.: Has primary production declined in the Salish Sea?, Can. J. Fish. Aquat. Sci., 363–6310, https://doi.org/10.1139/cjfas-2020-0115, 2020.
Jones, A. R., Hosegood, P., Wynn, R. B., De Boer, M. N., Butler-Cowdry, S., and Embling, C. B.: Fine-scale hydrodynamics influence the spatio-temporal distribution of harbour porpoises at a coastal hotspot, Prog. Oceanogr., 128, 30–48, https://doi.org/10.1016/j.pocean.2014.08.002, 2014.
Kärnä, T., Ljungemyr, P., Falahat, S., Ringgaard, I., Axell, L., Ko- rabel, V., Murawski, J., Maljutenko, I., Lindenthal, A., Jandt- Scheelke, S., Verjovkina, S., Lorkowski, I., Lagemaa, P., She, J., Tuomi, L., Nord, A., and Huess, V.: Nemo-Nordic 2.0: operational marine forecast model for the Baltic Sea, Geosci. Model Dev., 14, 5731–5749, https://doi.org/10.5194/gmd-14- 5731-2021, 2021.
Kendall, M. G.: Rank correlation methods, 1st edn., Griffin, London, vii + 196 p., https://webcat.library.ubc.ca/vwebv/ holdingsInfo?bibId=2313941 (last access: 14 January 2025), 1948.
Khangaonkar, T., Sackmann, B., Long, W., Mohamedali, T., and Roberts, M.: Simulation of annual biogeochemical cycles of nutrient balance, phytoplankton bloom(s), and DO in Puget Sound using an unstructured grid model, Ocean Dynam., 62, 1353– 1379, https://doi.org/10.1007/S10236-012-0562-4, 2012.
Khangaonkar, T., Nugraha, A., Xu, W., and Balaguru, K.: Salish Sea Response to Global Climate Change, Sea Level Rise, and Future Nutrient Loads, J. Geophys. Res.-Oceans, 124, 3876–3904, https://doi.org/10.1029/2018JC014670, 2019.
Khangaonkar, T., Nugraha, A., Lakshitha, P., Keister, J. E., and Borde, A.: Projections of algae, eelgrass, and zooplankton ecological interactions in the inner Salish Sea – for future climate and altered oceanic states, Ecol. Model., 441, 109420, https://doi.org/10.1016/j.ecolmodel.2020.109420, 2021a.
Khangaonkar, T., Nugraha, A., Yun, S. K., Premathilake, L., Keister, J. E., and Bos, J.: Propagation of the 2014–2016 Northeast Pacific Marine Heatwave Through the Salish Sea, Front. Mar. Sci., 8,	787604, https://doi.org/10.3389/fmars.2021.787604, 2021b.
Levier, B., Tréguier, A.-M., Madec, G., and Garnier, V.: Free surface and variable volume in the NEMO code, IFremer (Project Deliverable Report No. MERSEA-WP09-CNRS-STR-03-1A; p. 46, France, Zenodo, https://doi.org/10.5281/zenodo.3244182, 2007.
Libralato, S. and Solidoro, C.: Bridging biogeochemical and food web models for an End-to-End representation of marine ecosystem dynamics: The Venice lagoon case study, Ecol. Model., 220, 2960–2971, 2009.
Lin, Y., Dunphy, M., Krassovski, M., Blanken, H., and Hourston, R.: Preliminary development of a relocatable ocean model system for the west coast of Canada: A summary of OPP FA5 modeling (Canadian Technical Report of Hydrography and Ocean Sciences No. 343; p. v + 25 p.), Fisheries & Oceans Canada, ISBN: 978-0-660-45245-6, https://waves-vagues.dfo-mpo.gc. ca/library-bibliotheque/41073794.pdf (last access: 10 January 2025), 2022.
Loarie, S. R., Duffy, P. B., Hamilton, H., Asner, G. P., Field, C. B., and Ackerly, D. D.: The velocity of climate change, Nature, 462, 1052–1055, https://doi.org/10.1038/nature08649, 2009.
MacCready, P., McCabe, R. M., Siedlecki, S. A., Lorenz, M., Giddings, S. N., Bos, J., Albertson, S., Banas, N. S., and Garnier, S.: Estuarine Circulation, Mixing, and Residence Times in the Salish Sea, J. Geophys. Res.-Oceans, 126, e2020JC016738, https://doi.org/10.1029/2020JC016738, 2021.
Macias, D., Garcia-Gorriz, E., Piroddi, C., and Stips, A.: Biogeochemical control of marine productivity in the Mediterranean Sea during the last 50 years, Glob. Biogeochem. Cycles, 28, 897– 907, 2014.
Madec, G., Bourdallé-Badie, R., Bouttier, P.-A., Bricaud, C., Bru- ciaferri, D., Calvert, D., Chanut, J., Clementi, E., Coward, A., Delrosso, D., Ethé, C., Flavoni, S., Graham, T., Harle, J., Iovino, D., Lea, D., Lévy, C., Lovato, T., Martin, N., Masson, S., Mocavero, S., Paul, J., Rousset, C., Storkey, D., Storto, A., and Vancoppenolle, M.: NEMO ocean engine, Zenodo [code], https://doi.org/10.5281/zenodo.1464816, 2017.
Madec, G., Bell, M., Blaker, A., Bricaud, C., Bruciaferri, D., Cas- trillo, M., Calvert, D., Chanut, J., Clementi, E., Coward, A., Epicoco, I., Éthé, C., Ganderton, J., Harle, J., Hutchinson, K., Iovino, D., Lea, D., Lovato, T., Martin, M., Martin, N., Mele, F., Martins, D., Masson, S., Mathiot, P., Mele, F., Mocavero, S., Müller, S., Nurser, A. J. G., Paronuzzi, S., Peltier, M., Person, R., Rousset, C., Rynders, S., Samson, G., Téchené, S., Vancoppenolle, M., and Wilson, C.: NEMO Ocean Engine Reference Manual, Zen- odo, https://doi.org/10.5281/zenodo.8167700, 2023.
Mann, H. B.: Nonparametric Tests Against Trend, Econometrica, 13, 245–259, https://doi.org/10.2307/1907187, 1945.
Mantua, N. J. and Hare, S. R.: The Pacific Decadal Oscillation, J. Oceanogr., 58,	35–44, https://doi.org/10.1023/A:1015820616384, 2002.
Mantua, N. J., Hare, S. R., Zhang, Y., Wallace, J. M., and Francis, R. C.: A Pacific Interdecadal Climate Oscillation with Impacts on Salmon Production, B. Am. Me- teorol. Soc., 78, 1069–1079, https://doi.org/10.1175/1520- 0477(1997)078<1069:APICOW>2.0.CO;2, 1997.
Martins, E. G., Hinch, S. G., Patterson, D. A., Hague, M. J., Cooke, S. J., Miller, K. M., Lapointe, M. F., English, K. K., and Farrell, A. P.: Effects of river temperature and climate warming on stock-specific survival of adult migrating Fraser River sockeye salmon (Oncorhynchus nerka), Glob. Change Biol., 17, 99–114, https://doi.org/10.1111/j.1365-2486.2010.02241.x, 2011.
Masson, D.: Deep Water Renewal in the	Strait of Georgia, Estuar.	Coast.	Shelf Sci.,	54,	115–126, https://doi.org/10.1006/ecss.2001.0833, 2002.
Masson, D. and Cummins, P. F.: Temperature trends and interannual variability in the Strait of Georgia, British Columbia, Cont. Shelf Res., 27, 634–649, 2007.
Masson, D. and Peña, A.: Chlorophyll distribution in a temperate estuary: The Strait of Georgia and Juan de Fuca Strait, Estuar. Coast. Shelf Sci., 82, 19–28, https://doi.org/10.1016/j.ecss.2008.12.022, 2009.
McPhaden, M. J., Zebiak, S. E., and Glantz, M. H.: ENSO as an Integrating Concept in Earth Science, Science, 314, 1740–1745, https://doi.org/10.1126/science.1132588, 2006.
Melnychuk, M. C., Welch, D. W., and Walters, C. J.: Spatiotemporal migration patterns of Pacific Salmon smolts in rivers and coastal marine waters, PLoS ONE, 5, e12916, https://doi.org/10.1371/journal.pone.0012916, 2010.
Melnychuk, M. C., Christensen, V., and Walters, C. J.: Meso-scale movement and mortality patterns of juvenile coho salmon and steelhead trout migrating through a coastal fjord, Environ. Biol. Fishes, 96, 325–339, https://doi.org/10.1007/s10641-012-9976- 6, 2013.
Miller, K. M., Teffer, A., Tucker, S., Li, S., Schulze, A. D., Trudel, M., Juanes, F., Tabata, A., Kaukinen, K. H., Ginther, N. G., Ming, T. J., Cooke, S. J., Hipfner, J. M., Patterson, D. A., and Hinch, S. G.: Infectious disease, shifting climates, and opportunistic predators: Cumulative factors potentially impacting wild salmon declines, Evol. Appl., 7, 812–855, https://doi.org/10.1111/eva.12164, 2014.
Millero, F.: History of the Equation of State of Seawater, Oceanography, 23, 18–33, https://doi.org/10.5670/oceanog.2010.21, 2010.
Moore, S. K., Johnstone, J. A., Banas, N. S., and Salathé, E. P.: Present-day and future climate pathways affecting Alexandrium blooms in Puget Sound, WA, USA, Harmful Algae, 48, 1–11, https://doi.org/10.1016/j.hal.2015.06.008, 2015.
Morrison, J., Quick, M. C., and Foreman, M. G. G.: Climate change in the Fraser River watershed: flow and temperature projections, J. Hydrol., 263, 230–244, 2002.
Morrison, J., Foreman, M. G. G., and Masson, D.: A method for estimating monthly freshwater discharge affecting British Columbia coastal waters, Atmos. Ocean, 50, 1–8, https://doi.org/10.1080/07055900.2011.637667, 2012.
NOAA National Data Buoy Center: Meteorological and oceanographic data collected from the National Data Buoy Center Coastal-Marine Automated Network (C-MAN) and moored (weather) buoys, NOAA National Centers for Environmental Information [data set], https://www.ncei.noaa.gov/archive/ accession/NDBC-CMANWx (last access: 10 January 2025), 1971.
Oldford, G.: goldford/HOTSSea_v1_NEMOandSupportCode: HOTSSea v1 Code for GMD-2024-58 v2, Zenodo [code], https://doi.org/10.5281/zenodo.13887813, 2024a.
Oldford, G.: HOTSSea Model v1 monthly outputs (v1.02.4), Zen- odo [data set], https://doi.org/10.5281/zenodo.13942109, 2024b.
Oldford, G. and Dunphy, M.: HOTSSea Forcings – HRDPS, Zen- odo [data set] https://doi.org/10.5281/zenodo.12193924, 2024a.
Oldford, G. and Dunphy, M.: HOTSSea v1 Forcings – CIOPSW, ERA5, ORAS5, Obs, Runoff etc, Zenodo [data set], https://doi.org/10.5281/zenodo.12312769, 2024b.
Olson, E. M., Allen, S. E., Do, V., Dunphy, M., and Ian- son, D.: Assessment of Nutrient Supply by a Tidal Jet in the Northern Strait of Georgia Based on a Biogeochemical Model, J. Geophys. Res.-Oceans, 125, e2019JC015766, https://doi.org/10.1029/2019jc015766, 2020.
Pacific Salmon Foundation (PSF): Digital elevation model for the Salish Sea at a resolution of 80 metres, Pacific Salmon Foundation (PSF) [data set], https://doi.org/10.48689/6b05c683-10f5- 44fe-9983-7914e5cf0c72, 2014.
Pacific Salmon Foundation (PSF): CSOP- CTD readings in the Strait of Georgia [data set], v31-12-2022, https://soggy2.zoology.ubc.ca/geonetwork/srv/eng/catalog. search#/metadata/21cf5dad-692e-42ac-98e6-e9067289814c (last access: 14 January 2025), 2023.
Parzen, E.: An approach to empirical time series analysis, Radio Sci., 68, 937–951, 1964.
Pata, P. R., Galbraith, M., Young, K., Margolin, A. R., Perry, R. I., and Hunt, B. P. V.: Persistent zooplankton bioregions reflect long-term consistency of community composition and oceanographic drivers in the NE Pacific, Prog. Oceanogr., 206, 102849, https://doi.org/10.1016/j.pocean.2022.102849, 2022.
Pauly, D.: The gill-oxygen limitation theory (GOLT) and its critics, Sci. Adv., 7, eabc6050, https://doi.org/10.1126/sciadv.abc6050, 2021.
Pauly, D. and Cheung, W. W. L.: Sound physiological knowledge and principles in modeling shrinking of fishes under climate change, Glob. Change Biol., 24, e15–e26, https://doi.org/10.1111/gcb.13831, 2018.
Pawlowicz, R., Hannah, C., and Rosenberger, A.: Lagrangian observations of estuarine residence times, dispersion, and trapping in the Salish Sea, Estuar. Coast. Shelf Sci., 225, 106246, https://doi.org/10.1016/j.ecss.2019.106246, 2019.
Pearsall, A. I., Schmidt, M., Kemp, I., and Riddell, B.: Factors Limiting Survival of Juvenile Chinook Salmon, Coho Salmon and Steelhead in the Salish Sea: Synthesis of Findings of the Salish Sea Marine Survival Project, Pacific Salmon Fountain, Vancouver, BC, https://marinesurvival.wpengine.com/ wp- content/uploads/2021PSF- SynthesisPaper- Screen.pdf (last access: 14 January 2025), 2021.
Peña, M. A., Masson, D., and Callendar, W.: Annual plankton dynamics in a coupled physical–biological model of the Strait of Georgia, British Columbia, Prog. Oceanogr., 146, 58–74, 2016.
Perry, I.: Vignette 2: Lower Trophic Levels in the Salish Sea, in: State Salish Sea, edited by: Sobocinski, K. L., Broadhurst, G., and Baloy, N. J. K., Salish Sea Institute, Western Washington University, https://doi.org/10.25710/vfhb-3a69, 2021.
Piroddi, C., Akoglu, E., Andonegi, E., Bentley, J. W., Celiæ, I., Coll, M., Dimarchopoulou, D., Friedland, R., de Mutsert, K., Girardin, R., Garcia-Gorriz, E., Grizzetti, B., Hernvann, P. Y., Heymans, J. J., Müller-Karulis, B., Libralato, S., Ly- nam, C. P., Macias, D., Miladinova, S., Moullec, F., Palialexis, A., Parn, O., Serpetti, N., Solidoro, C., Steenbeek, J., Stips, A., Tomczak, M. T., Travers-Trolet, M., and Tsikliras, A. C.: Effects of Nutrient Management Scenarios on Marine Food Webs: A Pan-European Assessment in Support of the Marine Strategy Framework Directive, Front. Mar. Sci., 8, 596797, https://doi.org/10.3389/fmars.2021.596797, 2021.
Potisomporn, P., Adcock, T. A. A., and Vogel, C. R.: Evaluating ERA5 reanalysis predictions of low wind speed events around the UK, Energy Rep., 10, 4781–4790, https://doi.org/10.1016/j.egyr.2023.11.035, 2023.
Preikshot, D. B.: The influence of geographic scale, climate and trophic dynamics upon North Pacific oceanic ecosystem models, University of British Columbia, 208 pp., https://doi.org/10.14288/1.0074902, 2007.
Riche, O., Johannessen, S. C., and Macdonald, R. W.: Why timing matters in a coastal sea: Trends, variability and tipping points in the Strait of Georgia, Canada, J. Mar. Syst., 131, 36–53, 2014.
Riddell, B., Brodeur, R. D., Bugaev, A. V., Moran, P., Murphy, J. M., Orsi, J. A., Trudel, M., Weitkamp, L. A., Wells, B. K., and Wertheimer, A. C.: Ocean ecology of Chinook salmon, in: Ocean Ecology of Pacific Salmon and Trout, edited by: Beamish, R. J., American Fisheries Society, Bethesda, Maryland, 555–696, ISBN 978-1-934874-45-5, 2018.
Rose, K. A.: End-to-end models for marine ecosystems: Are we on the precipice of a significant advance or just putting lipstick on a pig?, Sci. Mar., 76, 195–201, https://doi.org/10.3989/scimar.03574.20B, 2012.
Rose, K. A., Allen, J. I., Artioli, Y., Barange, M., Blackford, J., Carlotti, F., Creekmore, S., Cropp, R., Daewel, U., Edwards, K., Flynn, K., Hill, S. L., HilleRisLambers, R., Huse, G., Mackinson, S., Mergrey, B., Moll, A., Rivkin, R., Salihoglu, B., Schrum, C., Shannon, L., Shin, Y., Smith, S. L., Smith, C., Solidoro, C., St. John, M., and Zhou, M.: End-To-End Models for the Analysis of Marine Ecosystems: Challenges, Issues, and Next Steps, Mar. Coast. Fish. Dyn. Manag. Ecosyst. Sci., 2, 115–130, https://doi.org/10.1577/C09-059.1, 2010.
Schulzweida, U.:	CDO User Guide, Zenodo, https://doi.org/10.5281/zenodo.7112925, 2022.
Seabold, S. and Perktold, J.: Statsmodels: Econometric and Statistical Modeling with Python, in: Proceedings of the 9th Python in Science Conference, Austin, Texas, 28 June–3 July 2010, 92–96, https://doi.org/10.25080/Majora-92bf1922-011, 2010.
Sen, P. K.: Estimates of the Regression Coefficient Based on Kendall’s Tau, J. Am. Stat. Assoc., 63, 1379–1389, https://doi.org/10.1080/01621459.1968.10480934, 1968.
Shapiro, S. S. and Wilk, M. B.: An Analysis of Variance Test for Normality (Complete Samples), Biometrika, 52, 591–611, https://doi.org/10.2307/2333709, 1965.
Sharma, R., Vélez-Espino, L. A., Wertheimer, A. C., Mantua, N., and Francis, R. C.: Relating spatial and temporal scales of climate and ocean variability to survival of Pacific Northwest Chinook salmon (Oncorhynchus tshawytscha), Fish. Oceanogr., 22, 14–31, https://doi.org/10.1111/fog.12001, 2013.
Snauffer, E. L., Masson, D., and Allen, S. E.: Modelling the dispersal of herring and hake larvae in the Strait of Georgia for the period 2007–2009, Fish. Oceanogr., 23, 375–388, https://doi.org/10.1111/fog.12072, 2014.
Sobocinski, K. L., Kendall, N. W., Greene, C. M., and Schmidt, M. W.: Ecosystem indicators of marine survival in Puget Sound steelhead trout, Prog. Oceanogr., 188, 102419, https://doi.org/10.1016/j.pocean.2020.102419, 2020.
Sobocinski, K. L., Greene, C. M., Anderson, J. H., Kendall, N. W., Schmidt, M. W., Zimmerman, M. S., Kemp, I. M., Kim, S., and Ruff, C. P.: A hypothesis-driven statistical approach for identifying ecosystem indicators of coho and Chinook salmon marine survival, Ecol. Indic., 124, 107403, https://doi.org/10.1016/j.ecolind.2021.107403, 2021.
Soontiens, N. and Allen, S. E.: Modelling sensitivities to mixing and advection in a sill-basin estuarine system, Ocean Model., 112, 17–32, https://doi.org/10.1016/j.ocemod.2017.02.008, 2017.
Soontiens, N., Allen, S. E., Latornell, D., Souëf, K. L., Paquin, J., Lu, Y., Thompson, K., Korabel, V., Soontiens, N., Allen, S. E., Latornell, D., Souëf, K. L., Paquin, J., Lu, Y., Thompson, K., Korabel, V., Surges, S., Soontiens, N., Allen, S. E., Latornell, D., Souëf, K. L., and Machuca, I.: Storm Surges in the Strait of Georgia Simulated with a Regional Model, Atmos. Ocean, 54, 1–21, https://doi.org/10.1080/07055900.2015.1108899, 2016.
Suchy, K. D., Young, K., Galbraith, M., Perry, R. I., and Costa, M.: Match/Mismatch Between Phytoplankton and Crustacean Zooplankton Phenology in the Strait of Georgia, Canada, Front. Mar. Sci., 9, 832684, https://doi.org/10.3389/fmars.2022.832684, 2022.
Teffer, A. K., Bass, A. L., Miller, K. M., Patterson, D. A., Juanes, F., and Hinch, S. G.: Infections, fisheries capture, temperature, and host responses: multistressor influences on survival and behaviour of adult Chinook salmon, Ocean Track. Netw. Adv. Aquat. Res. Manag., 75, 2069–2083, https://doi.org/10.1139/cjfas-2017-0491, 2018.
Theil, H.: A rank-invariant method of linear and polynominal regression analysis, Indagationes Mathematicae, 12, 1397–1412, 1950.
Thomson, R. E. and Huggett, W. S.: M2 Baroclinic Tides in Johnstone Strait, British Columbia, J. Phys. Oceanogr., 10,	1509–1539,	https://doi.org/10.1175/1520- 0485(1980)010<1509:MBTIJS>2.0.CO;2, 1980.
Tietsche, S., Alonso-Balmaseda, M., Zuo, H., and de Rosnay, P.: Comparing Arctic winter sea-ice thickness from SMOS and ORAS5, Tech. Memo., https://www.ecmwf.int/en/elibrary/ 17275- comparing- arctic- winter- sea- ice- thickness- smos- and- oras5 (last access: 14 January 2025), 2017.
Trudel, M. and Tucker, S.: Depth distribution of 1SW Chinook salmon in Quatsino Sound, British Columbia, during winter, North Pacific Anadromous Fish Commission, Research Document No. 1453, 8 pp., https://www.npafc.org/wp-content/ (last access: 10 January 2025), 2013.
Trudel, M., Fisher, J., Orsi, J. A., Morris, J. F. T., Thiess, M. E., Sweeting, R. M., Hinton, S., Fergusson, E. A., and Welch, D. W.: Distribution and Migration of Juvenile Chinook Salmon Derived from Coded Wire Tag Recoveries along the Continental Shelf of Western North America, Trans. Am. Fish. Soc., 138, 1369–1391, https://doi.org/10.1577/t08-181.1, 2009.
Tucker, S., Trudel, M., Welch, D. W., Candy, J. R., Morris, J. F. T., Thiess, M. E., Wallace, C., Teel, D. J., Crawford, W., Farley, E. V., and Beacham, T. D.: Seasonal Stock-Specific Migrations of Juvenile Sockeye Salmon along the West Coast of North America: Implications for Growth, Trans. Am. Fish. Soc., 138, 1458– 1480, https://doi.org/10.1577/T08-211.1, 2009.
Tuller, S. E.: Measured wind speed trends on the west coast of Canada, Int. J. Climatol., 24, 1359–1374, https://doi.org/10.1002/joc.1073, 2004.
Virtanen, P., Gommers, R., Oliphant, T. E., et al.: SciPy 1.0: fundamental algorithms for scientific computing in Python, Nat. Methods, 17, 261–272, https://doi.org/10.1038/s41592-019-0686-2, 2020.
Walters, C. J. and Christensen, V.: Effect of nonadditivity in mortality rates on predictions of potential yield of forage fishes, Ecol. Model., 410, 108776, https://doi.org/10.1016/j.ecolmodel.2019.108776, 2019.
Walters, C. J., Pauly, D., and Christensen, V.: Ecospace: Prediction of Mesoscale Spatial Patterns in Trophic Relationships of Exploited Ecosystems, with Emphasis on the Impacts of Marine Protected Areas, Ecosystems, 2, 539–554, https://doi.org/10.1007/s12009-000-0069-3, 1999.
Walther, B. A. and Moore, J. L.: The concepts of bias, precision and accuracy, and their use in testing the performance of species richness estimators, with a literature review of estimator performance, Ecography, 28, 815–829, https://doi.org/10.1111/j.2005.0906-7590.04112.x, 2005.
Welch, D. W., Melnychuk, M. C., Payne, J. C., Rechisky, E. L., Porter, A. D., Jackson, G. D., Ward, B. R., Vincent, S. P., Wood, C. C., and Semmens, J.: In situ measurement of coastal ocean movements and survival of juvenile Pacific salmon, P. Natl. Acad. Sci. USA, 108, 8708–8713, https://doi.org/10.1073/pnas.1014044108, 2011.
White, H.: A Heteroskedasticity-Consistent Covariance Matrix Estimator and a Direct Test for Heteroskedasticity, Econometrica, 48, 817–838, https://doi.org/10.2307/1912934, 1980.
Wilcox, R.: A Note on the Theil-Sen Regression Estimator When the Regressor Is Random and the Error Term Is Heteroscedas- tic, Biom. J., 40, 261–268, https://doi.org/10.1002/(SICI)1521- 4036(199807)40:3<261::AID-BIMJ261>3.0.CO;2-V, 1998.
Willmott, C. J.: On the Validation of Models, Phys. Geogr., 2, 184– 194, https://doi.org/10.1080/02723646.1981.10642213, 1981.
Yin, K., Goldblatt, R. H., Harrison, P. J., John, M. A. S., Clifford, P. J., Beamish, R. J., St. John, M. A., Clifford, P. J., and Beamish, R. J.: Importance of wind and river discharge in influencing nutrient dynamics and phytoplankton production in summer in the central Strait of Georgia, Mar. Ecol. Prog. Ser., 161, 173–183, https://doi.org/10.3354/meps161173, 1997.
Yue, S., Pilon, P., Phinney, B., and Cavadias, G.: The influence of autocorrelation on the ability to detect trend in hydrological series, Hydrol. Process., 16, 1807–1829, https://doi.org/10.1002/hyp.1095, 2002.
Zaron, E. D.: Introduction to Ocean Data Assimilation, in: Operational Oceanography in the 21st Century, edited by: Schiller, A. and Brassington, G. B., Springer Netherlands, Dordrecht, 321– 350, https://doi.org/10.1007/978-94-007-0332-2_13, 2011.
Zuo, H., Balmaseda, M. A., Tietsche, S., Mogensen, K., and Mayer, M.: The ECMWF operational ensemble reanalysis–analysis system for ocean and sea ice: a description of the system and assessment, Ocean Sci., 15, 779–808, https://doi.org/10.5194/os- 15-779-2019, 2019.
The Salish Sea is an inland sea in the northeast Pacific spanning Canadian and American waters with estuarine characteristics, fjords, and high marine biodiversity (Harrison et