An overview of the NEMO modelling for the BaySys project

The Hudson Bay Complex (HBC) is a large shallow inland sea in north central Canada. It includes Hudson Bay (HB) itself, as well as the neighbouring James Bay, Foxe Basin, Hudson Strait and Ungava Bay. The main connection to the global ocean is to the North Atlantic through Hudson Strait, but there is also a connection to the Canadian Arctic through the small and narrow Fury and Hecla Strait between Foxe Basin and the Gulf of Boothia. The region is ice-covered for much of the year, although the time period of ice coverage has been decreasing in recent decades (Hochheim and Barber, 2010; Hochheim et al., 2011; Castro de la Guardia et al., 2017). The basin receives freshwater from precipitation (McCullough et al., 2019), as well as from river runoff from a large watershed that provides over 900 km³ of freshwater per year based on a draining area of 3.8 × 10⁶ km² (Shiklomanov and Shiklomanov, 2003; Ridenour et al., 2019b). Thus, the bay is quite fresh and generally well stratified (Ingram and Prinsenberg, 1998). The presence of sea ice and the large volume of freshwater creates a critical habitat for a large number of species (Stewart and Lockhart, 2004). Despite this importance, the bay has not been studied in as much detail as other water bodies across Canada. Changes in continental flow regimes have also been imposed by increasing anthropogenic development and river control due to regulation (Dynesius and Nilsson, 1994), defined here as the control of river stage or discharge for human consumption or energy production needs. Greater than 47% of the freshwater entering the bay is being regulated intensely (>70% is considered regulated), which impacts the timing and magnitude of discharge. The BaySys project was created to provide a scientific basis to separate climate change effects from those of regulation of freshwater on physical, biological and biogeochemical conditions in HB. The project addresses this objective from a “systems” perspective, with sub-objectives to examine the climate, marine and freshwater systems and to study the cycling of carbon and contaminants. This introduction article is intended to introduce and discuss the development of a modelling framework to examine simulated climate change and river discharge regulation and their combined impact on marine conditions in the HBC.

Regional general circulation models (GCMs) were first used in the early 1990s to study the sea ice (Wang et al., 1994a) and ocean (Wang et al., 1994b) circulation in HB. Another early modelling study (Saucier and Dionne, 1998) included the role of regulated runoff from hydroelectric development, albeit in a simple idealized sense, to examine the factors impacting the seasonal cycle of low salinity waters within the bay. Additionally, the runoff was applied in terms of a simple salt flux function rather than an actual freshwater flux (Saucier and Dionne, 1998). Meanwhile, Gough (1998) used a simple 1-D column model, as well as a coarse resolution global ocean model, to predict that sea level rise due to global warming may exceed the local sea level decrease from isostatic rebound within the next century. Output from climate models began to be applied to HBC questions during this time period (Ingram et al., 1996; Gagnon and Gough, 2005).

However, the coarse resolution used in the climate models of the time meant that forced ocean/sea ice models continued to be used to examine questions of climate sensitivity because of their higher resolution (Gough and Allakhverdova, 1999; Gough, 2001). Increased resolution and better representation of mixing, including tidal forcing, allowed Saucier et al. (2004) to reproduce many features of the seasonal circulation and hydrography in the HBC over a 2-year integration period. This same model was later used by St-Laurent et al. (2008) to examine the role of under-ice friction on the tides and by Defossez et al. (2012) to compare winter and summer circulation in Foxe Basin. St-Laurent et al. (2011) then showed that the annual freshwater balance is mainly between river runoff and export through Hudson Strait, albeit with substantial seasonal exchange of freshwater into the interior of the basin through Ekman transport. St-Laurent et al. (2012) also built a simple conceptual ocean model to further explore and understand issues related to freshwater content and export in the HBC. Joly et al. (2011) forced the model of Saucier et al. (2004) with projections of future climate from a third generation coupled climate model (albeit without considering changes in the hydrological forcing). Compared to a present day control run over 2001–2005, Joly et al. (2011) found that the greatest changes in sea ice and ocean heat content occur in James Bay, southeast Hudson Bay and Hudson Strait. Castro de la Guardia et al. (2013) used a sea ice model coupled to a slab ocean, driven by multiple 21st-century greenhouse gas emission scenarios, to find that spring sea ice concentration from 2090 to 2100 would be between 20% and 84% of present day. Recently, Kleptsova and Pietrzak (2018) used a coastal circulation model to show that tidal seasonality has changed significantly in the Canadian Arctic, including HB, over the past decades, due to the decaying extent of the sea ice. Meanwhile, Lemieux et al. (2018) included tides in a pan-Arctic coupled ice-ocean Nucleus for European Modelling of the Ocean (NEMO) simulation. Their domain of analysis was the Canadian Arctic Archipelago, including Foxe Basin, and they showed that tides were important in reducing excess landfast ice through increasing the ocean-ice stress.

Based on this previous work, the goal of the marine system modelling team within BaySys was to develop an integrated observational-modelling framework that can provide insights on, and improved representation of, physical, biological and biogeochemical processes in the HB system. The modelling is to provide a framework and tool for simulating projected changes in marine state and dynamic variables, while also enabling an integration of observations and numerical analyses. Given its widespread use by the Canadian research community, this modelling perspective is based on the NEMO ocean general circulation model coupled to version 2 of the Louvain-la-Neuve sea ice model (LIM2). Thus, we use an Arctic NEMO configuration (Figure 1a) to provide a local (20–100 km, estuarine and coastal) and a regional (100–1,000 km, bay-wide) perspective and understanding of changes in freshwater-marine coupling in response to a changing climate in both regulated and naturalized regimes. The Arctic configuration further enables a link between Arctic and sub-Arctic domains and ocean, sea ice, and atmospheric phenomena to consider the tightly integrated nature of the high-latitude climate system in examining impacts on the HBC. We note that the pan-Arctic domain is essential in ensuring that the climate change signal (a hemispheric-scale phenomenon) within the HBC is simulated adequately. Previous studies (Ingram and Prinsenberg, 1998; Jones et al., 2003) have demonstrated that waters from the Canadian Arctic, Siberian rivers and Pacific Water (via Bering Strait) all enter the HBC over timescales of 2 years to 10 years. Given that HB is filled with Arctic waters and that climate change is expected to have the largest impacts at high latitudes, understanding the response of the HBC to a changing climate requires inputs that represent how the Arctic is responding to the climate forcing, potentially modifying the exchange into the HBC. Given the timescale of the response, the overall change to forcing and runoff from the Arctic will be the important changes and not short-term changes in river regulation in that region, as those effects will be integrated (minimized) over the transit times of the given waters to the HBC.

Work with existing NEMO configurations and results from experiments carried out early during the BaySys project are described in Section 2. Details of the main NEMO model configuration used, with the developments needed for the long climate experiments of the BaySys project, as well as a summary of the model forcing, are given in Section 3. Section 3 also includes a discussion and summary of the experiments carried out and their setup. Section 4 provides some evaluation of the BaySys NEMO experiments with the BaySys moorings, as well as satellite data. Section 5 finishes with a discussion of the development of the BaySys NEMO configuration, including limitations, as well as thoughts on where the research could be taken in the future.

While the main NEMO model configuration, described in Section 3, was being developed and implemented for the long climate experiments, existing NEMO configurations were used to carry out sensitivity experiments and begin to understand aspects of the circulation of the HBC. Some of the setups were carried out within the framework of BaySys, while other studies were not, but all helped to provide insight into the model and its functioning in the HBC.

Ridenour et al. (2019b) considered present-day freshwater dynamics in the HBC, in addition to evaluating the sensitivity of the region to model resolution and runoff forcing. Using different estimates of runoff allowed this work to analyse the model sensitivity to the amount, seasonality and timescale of the runoff product. While surface fluxes, which included precipitation and evaporation as well as sea ice growth and melt, dominated the freshwater budget during the year, their net effect was small. The effect of river discharge, however, was large on an annual timescale and was offset by the nearly equal amount of freshwater advected out of the region. Decreased (increased) discharge from the data sets used resulted in weaker (stronger) circulation patterns within HB, which led to lower (higher) volume and freshwater exchange between the North Atlantic and HBC. Using higher model resolution, which also impacted the dynamics in HB, highlighted one particular region, between Southampton and Baffin Islands, where small-scale processes impacted exchange between Foxe Basin and Hudson Strait.

Historically, the strong geostrophic, cyclonic flow was a defining feature of HB that was found and supported by numerous observational and modelling studies. However, other studies (Gough et al., 2005; St-Laurent et al., 2011) hinted that this circulation may not be as stable as previously thought. Using a high-resolution ocean general circulation model, Ridenour et al. (2019a) showed the presence of a weak anticyclonic flow in eastern HB in summer. This finding was supported by satellite-based observations of absolute dynamic topography and geostrophic velocity. This flow, while geostrophic, is strongest through the centre of the bay and is induced by the spring freshet and strengthened by anticyclonic seasonal wind patterns.

The velocity, temperature and salinity fields from the Ridenour et al. (2019b) study were used to run a particle tracking model designed for ichthyoplankton (fish planktonic life stage; Schembri et al., 2021). Using particle backtracking in conjunction with otolith (fish ear bone) analysis, modelling fields were able to provide possible spawning temperatures and salinities of arctic cod in the HBC, aiding in proving the use of under-ice freshwater refugia in the Arctic for arctic cod survival (Schembri et al., 2021).

Jafarikhasragh et al. (2019) provided an assessment of the sensitivity of simulated sea surface temperature (SST) to atmospheric forcing and model resolution in the HBC. This study led to an improved understanding of bulk heat flux parameterizations in the NEMO model and how the model produces heat fluxes and drives SST on a basin-wide scale, with implications for air-sea heat flux characterization. Investigation of simulated thermodynamic and dynamic contributions to changes in sea ice thickness on seasonal timescales showed that thermodynamic processes govern sea ice thickness changes in the HBC (Jafarikhasragh, 2019). It further demonstrated that surface energy rather than ocean heat flux contributes to thermodynamic changes, while wind stress is associated with dynamic sea ice thickness changes. In light of demonstrated correspondence between observed and simulated sea ice conditions, results from this analysis led to an improved understanding of processes governing changes in sea ice thickness on seasonal timescales in the HBC, with implications for prediction.

Eastwood et al. (2020) examined processes controlling winter stratification in southeast HB. Although mainly an observation-based study, they (Eastwood et al., 2020) used NEMO model output to help build a conceptual model that had James Bay river inflow driving northward transport of both river water and brines past the Belcher Islands.

Tao and Myers (2021) used a high-resolution NEMO configuration to highlight circulation pathways of potential pollutants due to ocean advection, through the use of the Ariane Lagrangian particle tracking tool (Blanke and Raynaud, 1997; Blanke et al., 1999). They (Tao and Myers, 2021) found three main pathways of circulation for particles released in the HBC: (1) eastward outflow via the southern coast of Hudson Strait into the Labrador Current; (2) inflow with the Baffin Island Current along the northern coast of Hudson Strait; and (3) the cyclonic circulation in Foxe Basin and HB. This last point implied that for particles released in the interior of the HBC, they mainly remained in the HBC.

Castro de la Guardia et al. (2019) used NEMO and the biogeochemistry (BGC) with Light Iron and Nutrient limitation and Gases (BLING) to examine the linkages between sea ice loss and high-frequency winds on Arctic net primary production and biogenic carbon export. They (Castro de la Guardia et al., 2019) examined a number of regions around the pan-Arctic, but the HBC showed that high-frequency winds are crucial for generating the autumn bloom in seasonally ice-covered seas. The high-frequency winds were found to be responsible for over 50% of the autumn primary production and carbon export and between 20% and 40% of the same fields in spring. This same work also included a detailed model evaluation relevant to their biogeochemical modelling study and showed, in the HBC, very good correspondence of the simulated NEMO sea ice fields, averaged over the bay, to passive microwave data (Comiso, 2000). They also showed that the model-simulated chlorophyll-a (chl-a) seasonality agreed well with ocean colour measurements from the Moderate Resolution Imaging Spectroradiometer sensors on board the Aqua satellite (MODIS-A; Ocean Biology Processing Group, 2003), albeit with slight overestimates of the magnitudes of the spring and autumn blooms. Finally, the observed estimate of early fall mean primary production is approximately 200 mg C m⁻² day⁻¹ (derived from 12 stations in figure 4a in Lapoussière et al., 2013), which is only 50 mg C m⁻² day⁻¹ lower than the simulated regional autumn mean.

The model used is the NEMO numerical framework version 3.6 (Madec et al., 2008). NEMO was chosen because it is widely used by the academic community, and is consistent with the framework developed by Canada’s intra-agency (Departments of Fisheries and Oceans, Environment and Climate Change and National Defence) Canadian Operational Network for Coupled Environmental Prediction Systems (CONCEPTS) initiative to allow direct integration of the knowledge gained in this project into ongoing governmental forecasting initiatives. The model is coupled with version 2 of the Louvain-la-Neuve sea ice model (LIM2; Fichefet and Morales Maqueda, 1997). With the need to consider the pan-Arctic domain, the configuration used is the Arctic and Northern Hemisphere Atlantic (ANHA), which uses a tri-polar grid with open boundaries at Bering Strait and 20°S in the Atlantic Ocean. These boundaries were chosen to be far enough from the study region to limit the impact of far-field behaviour on the HBC. Details on the data used at the boundaries, as well as sensitivity experiments to confirm their limited impact, are discussed in the following subsections. A nominal horizontal resolution of 0.25 degree (hereafter ANHA4; Figure 1a) is used to balance the need to represent boundary currents and some mesoscale processes while also allowing multiple century-long integrations to be carried out and passive tracers and biogeochemical model components to be included. The ANHA4 configuration has a resolution of 10–16 km in the HBC (Figure 1b). The configuration has 50 vertical levels, with the layer thickness smoothly increasing from 1.05 m at the surface to 453.13 m at the last level. A partial step is also enabled to better resolve the bathymetry. No temperature or salinity restoring is applied to avoid damping the runoff and climate signals we wish to study. Vertical mixing at subgrid scales is parameterized using the NEMO turbulent kinetic energy closure model (Madec et al., 2008). For lateral mixing, the model uses a bi-Laplacian operator with an eddy viscosity of 1.5 × 10¹¹ m⁴ s⁻¹. Subgrid-scale tracer lateral diffusion is parameterized with an isopycnal Laplacian operator with an iso-neutral eddy diffusivity of 300 m² s⁻¹, as typically used when running 1/4-degree NEMO experiments (e.g., Biastoch et al., 2021). The model baroclinic timestep is 1,080s. Tides are included by specifying geopotential tidal forcing with nine constituents in the momentum equations (K1, K2, M2, M4 N2, O1, P1, Q1, S2), as well as at the lateral open boundaries. The tidal constituents are taken from the TOPEX/Poseidon Global Inverse Solution TPXO 8 of the Oregon State University, USA (Egbert and Erofeeva, 2002). Using tides means the experiments are run with the nonlinear free surface and variable volume formulations of NEMO v3.6 activated. Initial conditions are from the Polar Hydrographic Climatology version 3.0 (PHC3.0; Steele et al., 2001).

A five-member ensemble of Coupled Model Intercomparison Project phase 5 (CMIP5) model experiments was chosen to provide the atmospheric fields to drive the ocean/sea ice model. The members were chosen to bracket future changes in temperature and precipitation across the HBC domain, thus maximizing climatic variability (Braun et al., 2021). The ensemble size was limited to five because of the expense of running 100-year-long experiments with NEMO. The choice of the CMIP5 model experiment was also impacted by the need to provide the required fields over the NEMO domain (Braun et al., 2021). Three general circulation models (GCM) were chosen. The Model for Interdisciplinary Research on Climate (MIROC5) was developed cooperatively by the Japanese research community for the IPCC Fifth Assessment (Watanabe et al., 2010) with resolution T85 for the atmosphere and the equivalent 1.4 degrees for the ocean except near the Equator (Watanabe et al., 2010). The Meteorological Research Institute Coupled General Circulation Model (MRI-CGCM3) was also developed in Japan for the Fifth Assessment, with atmospheric resolution of T159 and oceanic horizontal resolution of 1 degree in longitude and 0.5 in latitude (Yukimoto et al., 2012). The Geophysical Fluid Dynamics Laboratory Community Model (GFDL-CM3) was developed at the NOAA Geophysical Fluid Dynamics Laboratory with a cube-sphere grid, giving a C48 horizontal resolution (163 km to 231 km grid size; Donner et al., 2011). For both MIROC5 and MRI-CGCM3, two Representative Concentration Pathways (RCPs; IPCC, 2013), 4.5 and 8.5, were chosen. For the fifth member, an RCP4.5 experiment from GFDL-CM3 was used. Further details on the ensemble member selection are given in Stadnyk et al. (2019) and Braun et al. (2021).

To drive NEMO, the following fields were used from the GCM experiments: 3-hourly surface air temperature, atmospheric humidity, zonal and meridional surface winds and radiation (downwelling shortwave and longwave). Six-hourly surface atmospheric pressure was also needed for the biogeochemical modules discussed below. Temperature, precipitation and wind were bias-corrected using Watch Forcing Data, ERA-Interim (WFDEI), over the ocean domain, as discussed in Stadnyk et al. (2019) and Braun et al. (2021). Forcing fields were interpolated onto the NEMO model grid during the runs using the NEMO on-the-fly interpolation function (Madec et al., 2008). The Coordinated Ocean-ice Reference Experiments bulk formulae were applied to compute fluxes of heat, water, and momentum (Large and Yeager, 2009) at each model timestep. Monthly averaged boundary conditions at the model open boundaries of Bering Strait and 20°S were also taken from the output of the given CMIP5 experiment. Analysis of the baseline atmospheric conditions over the HBC for 2016–2018 is given in Lukovich et al. (2021c).

An additional historical control experiment was run from 1980 to 2018 using ERA-Interim forcing (Dee et al., 2011) and historical runoff, as discussed in the following subsection. Otherwise, this run was set up similarly to the historical phase of the climate model experiments. This experiment is defined as EPM017 in the figures and tables of this article.

Given the focus in BaySys on studying freshwater dynamics within the HBC, runoff forcing was an important aspect of the modelling experiments. Within the HBC, hydrological simulations were performed with a modified version of Arctic-HYdrological Predictions for the Environment (Arctic-HYPE; Andersson et al., 2015; Gelfan et al., 2017), a hydrological model that was improved and calibrated for the region (MacDonald et al., 2018; Stadnyk et al., 2020). Arctic-HYPE was forced by the same bias-corrected atmospheric forcing sets used to drive the NEMO simulations described above. Although MacDonald et al. (2018) used additional GCM simulations, they were not considered here because of the numerical expense of running century-long NEMO experiments. In all cases, monthly river discharge for the HBC was produced for each GCM/RCP pair, using the GCM historical forcing simulation for 1980–2005 and the future simulation for 2005–2070. Two configurations of Arctic-HYPE were run for each climate simulation, one naturalized scenario and one including river regulation (which reverses land cover changes by reservoir flooding and reverts lake outflow to pre-development conditions) in the HB domain. Thus, two sets of 90-year-long hydrological discharge scenarios were produced to drive NEMO-naturalized and regulated scenarios for each bias-corrected GCM/RCP pair. Additionally, historical WFDEI fields were used to produce naturalized and regulated runoff from 1980 to 2018 to drive a historical control simulation for comparison.

Given the dependence of the HBC on inputs from the Arctic, modelled hydrological scenarios were deemed necessary for the Arctic watersheds as well. Thus, an additional set of Arctic-HYPE simulations were carried out (Stadnyk et al., 2021) for the pan-Arctic domain, again driven by the same five GCM/RCP forcing sets for 1980–2070, plus the WFDEI historical forcing over 1980–2018. However, Arctic-HYPE was run with bias correction only for temperature and precipitation. Given a lack of details on regulation on Russian rivers, only naturalized output was produced for the pan-Arctic domain (Stadnyk et al., 2021).

HYPE model outputs satisfying the constraints of the BaySys project were not available for the rest of the model domain. Thus, river runoff was taken from the Canadian Centre for Climate Modelling and Analysis (CCCMA) CanESM2 model (Arora et al., 2011), based on historical (1950–2005), RCP4.5 (2006–2070) and RCP8.5 (2006–2070) experiments and a variable velocity flow river-routing algorithm (Arora and Boer, 1999).

Runoff was regridded from the river mouth positions onto the NEMO model grid using the approach discussed in Hu et al. (2018) and Hayashida et al. (2019). Within the same river mouth polygons, enhanced vertical mixing of 2 × 10⁻³ m² s⁻¹ was applied up to 30 m (i.e., constrained by water depth) of the water column to prevent unrealistically low salinities in long narrow estuaries, such as the Ob River.

Besides river runoff, freshwater is added into the high-latitude ocean by discharge from the Greenland Ice Sheet. This discharge has two components. Liquid melt, including tundra discharge, is added to the model similarly to river runoff. Solid discharge, or calving, is included in the model through a Lagrangian iceberg module (Marzocchi et al., 2015; Marson et al., 2018). This module includes the modification to apply ocean fields vertically through the thickness of the iceberg, as discussed in Marson et al. (2018).

For both liquid discharge and solid calving, we used fields from Lenaerts et al. (2015). As discussed in Gillard (2021), for the historical period, the Greenland Ice Sheet solid ice discharge was constructed from remote sensing records for 2000–2012 (Enderlin et al., 2014). A near-constant value based on this analysis period is then used to drive NEMO. Meanwhile, the liquid runoff portion of the Greenland Ice Sheet freshwater forcing originates from the runoff from Regional Atmospheric Climate Model version 2.1 (hereafter RACMO2.1; van Meijgaard et al., 2008). RACMO2.1 has a spatial resolution of approximately 11 km, is forced by ERA-Interim fields at its lateral boundaries, has a Greenland Ice Sheet surface mass balance (van Angelen et al., 2014), and improvements for the climate over Greenland (Ettema et al., 2010). Runoff is given spatial variability by the subdivision into eight basins. The historical scenario calculates runoff based on RACMO2.1 (1960–2012) for each basin. For the meltwater calculations beyond 2012, the regional climate model is forced with an atmospheric circulation climate model HadGEM2-ES. Runoff is distributed evenly to the ocean grid points along each individual basin and assimilated into the coupled land-atmosphere-ocean climate model Community Earth System Model version 1.1.2 (CESM1). The CESM1 is used to simulate multiple scenarios, three of which have been used in this study: a historical (1850–2005) and two future climate scenarios (2006–2200). Further details on how the RACMO fields are used in the NEMO model can be found in Gillard (2021). Gillard (2021) also showed that over the 2004–2016 period, the use of the Lenaerts et al. (2015) fields in NEMO leads to comparable results over the subpolar North Atlantic when compared to the more observationally based product of Bamber et al. (2012, 2018).

Previously, the ANHA4 configuration had been run coupled with the BLING version 0 (Galbraith et al., 2010), as discussed in Section 2. BLING version 0 is a reduced complexity biogeochemical model with four prognostic tracers: inorganic phosphate, dissolved organic phosphate, oxygen and iron. It diagnoses chl-a, phytoplankton production and particle export, considering light, macronutrient and iron limitations as well as temperature dependency. Using BLING as the choice for biogeochemical modelling has the benefit of lower computational demands, an advantage for running long simulations on high-resolution grids like ANHA4 (Castro de la Guardia et al., 2019). A detailed discussion of the setup and evaluation of ANHA4 and BLINGv0 can be found in Castro de la Guardia (2018) and Castro de la Guardia et al. (2019). Despite having only four prognostic tracers, BLING is able to reproduce the basic bloom dynamics and magnitude within the HBC (Castro de la Guardia, 2018; Castro de la Guardia et al., 2019). Given the interest in the carbon system within BaySys, and the need to use BLING while the BioGeoChemical Ice Incorporated Model (BIGCIIM) was being coupled to NEMO (see below for further details), BLING was updated to the newer BLINGv0 + DIC version (Galbraith et al., 2015), which adds dissolved carbon and alkalinity as prognostic variables (and a suite of diagnosed quantities).

For the purposes of BaySys, BLINGv0 + DIC was run coupled with NEMO for all three historical scenarios (MRI, MIROC, GFDL) over the historical period of 1980–2005 for both regulated and naturalized river runoff. For future periods, various practical considerations limited the future experiments with BLING to just the RCP8.5 scenarios of MRI and MIROC. In each case, the 2006–2070 period was run for each forcing for both naturalized and regulated runoff.

Biological fields were initialized with both observed climatology and model output. Dissolved oxygen and inorganic phosphate fields were derived from observed climatologies from World Ocean Atlas 2012 version 2 (WOA13; Garcia et al., 2014a, 2014b). Dissolved iron and organic phosphate come from the Geophysical Fluid Dynamics Laboratory (GFDL) Earth System Model version 2 (ESM2M) coupled with BLING (Galbraith et al., 2015). The GFDL ESM2M simulation is a global configuration at 1-degree nominal resolution and geopotential vertical coordinates. The simulation has a 100-year spin-up period using 1,860 forcing and an atmospheric carbon dioxide partial pressure (pCO₂) of 286 ppm. The initial conditions were built using the average of the last 20 years of the spin-up period.

The initial conditions of total alkalinity and dissolved inorganic carbon (DIC) derived from observed climatology are from the mapped product of GLobal Ocean Data Analysis Project version 2 (DICGLODAPv2; Key et al., 2015; Lauvset et al., 2016; Oslen et al., 2016). These fields were remapped onto the ANHA4 grid with units of mol m⁻³. However, DIC initial conditions (DIC_ic) were normalized to the simulation start year (DIC_ic = DIC_GLODAPv2 − DIC_diff). The variable DIC_diff is the anthropogenic carbon using estimates from DeVries (2014) and is calculated as the difference between DIC_DeVries(yri) and DIC_{DeVries(current – year)}.

Open boundary conditions for all tracers in BLING and atmospheric pCO₂ are derived from yearly output of the Community Earth System Model version 1.1.2 (CESM1) with biogeochemistry (BGC) simulations that were also part of the CMIP5 ensemble. We used the CESM1-BGC output 1 of ensemble r1i1p1 of the pre-industrial control and 20th-century experiments (RCP4.5 and RCP8.5; Lindsay et al., 2014). Tracer data were extracted at 20°S in the Atlantic Ocean and at Bering Strait (see Figure 1), the boundaries of the ANHA4 configuration. A source of iron at the ocean’s surface was added following the relation between dust deposition and iron concentrations described in Galbraith et al. (2010). The climatological monthly dust deposition input at the surface of the ANHA4 domain was derived from the Global Ozone Chemistry Aerosol Radiation and Transport model (Ginoux et al., 2001).

One limitation of BLING is that it only considers the pelagic system (plankton within the water column), yet the HBC experiences a seasonal ice cover. In terms of biogeochemistry, including/investigating the dynamics of both the sympagic (organisms associated with sea ice) and the pelagic systems and their interactions is important. Additionally, the limiting nutrient in the Arctic and sub-Arctic regions is nitrogen (Tremblay and Gagnon, 2009; Tremblay et al., 2009). Due to the importance of these factors, the plan was to couple the BioGeoChemical Ice Incorporated Model (BiGCIIM), based on the original model of Sibert et al. (2010, 2011), to NEMOv3.6 and LIM2 for BaySys. The prognostic tracers within BiGCIIM are primary producers (micro-algae), which are split into ice algae and a large and a small group of phytoplankton (diatoms and flagellates, respectively), secondary consumers (split into ice fauna, mesozooplankton and microzooplankton), particulate organic matter, dissolved organic matter, nitrate, ammonium and a bottom storage compartment. Mass balance is maintained through rates connecting each prognostic tracer.

A simple carbon module was then added to BiGCIIM (Lavoie et al., 2021). The carbon module calculates DIC and total alkalinity through the use of rates of respiration, production (uptake of carbon) and remineralization based on the original nitrogen model and the exchange between the ocean surface and atmosphere and DIC input from the rivers (Holliday et al., 2020).

The initial conditions for BiGCIIM were created from observational data available from the World Ocean Atlas 2018 (Garcia et al., 2018). Observational data were interpolated over the ANHA4 domain. These fields were used to create the open boundary condition for the model domain.

The river nutrient inputs are locally sensitive. Average nitrate, ammonium, particulate organic matter, dissolved organic matter, phosphate, dissolved inorganic carbon and total alkalinity were obtained from the observations collected by the Arctic Great Rivers Observatory (Holliday et al., 2020). To produce boundary conditions at the scale of the whole domain, the average nutrient inputs were then applied over the whole domain and multiplied by the discharge provided by the HYPE model regulated and naturalized runoff forcings, for both BLING and BiGCIIM. Future work will address nutrient input variability at the scale of individual hydrological basins.

Due to the time constraints linked to BaySys deadlines and the extent of work required to get BiGCIIM working, the decision was made to run BLINGv0-DIC with the physical oceanographic runs to have the availability of biogeochemical fields to analyse earlier and compare with BiGCIIM in the future. BiGCIIM will be run for all three historical scenarios (MRI-CGCM3, MIROC5, GFDL-CM3) over 1980–2005 for both regulated and naturalized river runoff. For the future periods, the RCP8.5 scenarios of MRI-CGCM3 and MIROC5 ideally will be run for 2006–2070 for each forcing for both naturalized and regulated runoff. Validation of BLINGv0+DIC for HBC is summarized in Section 4.

The BaySys NEMO experiments used to study the impacts of climate change and river regulation are summarized in the schematic in Figure 2. The historical period is defined as 1980–2005. Three pairs of experiments were run from 1980 to 2005, each pair having NEMO forced by the bias-corrected atmospheric forcing associated with the historical control run of each GCM (MIROC5, MRI-CGCM3 and GFDL-CM3). Each pair then includes A-HYPE river runoff produced using the forcing from the given bias-corrected historical GCM simulation. Additionally, a pair of historical control experiments were run using ERA-Interim forcing, with naturalized and regulated A-HYPE river runoff driven by ERA-based WFDEI forcing. The historical experiment with naturalized river runoff was then extended to 2018 to allow direct comparison with the BaySys observations and mooring records, as summarized in the following section.

Each historical experiment was then continued from the start of 2006 to 2070 using the climate forcing from each GCM. For MIROC5 and MRI-CGCM3, two RCP pairs (RCP4.5 and RCP8.5) were used. For GFDL-CM3, only RCP4.5 forcing was used. The end result is a 10-member ensemble of NEMO experiments over 2006–2070, five members each with naturalized and regulated river runoff forcing from A-HYPE.

To help evaluate NEMO, output from a historical control experiment was compared with different observational data sets. These include both BaySys moorings, sea ice data from various sources, and tide/water level data at Churchill.

Various fields, including temperature, salinity, currents and sea ice thickness, were taken from three long-term oceanographic moorings deployed as part of BaySys. These moorings were deployed in western (AN01, 59°58′N, 91°57′W; NE03, 57°50′N, 90°53′W) and southeastern (JB02, 54°41′N, 80°11′W) Hudson Bay (Figure 1b). The deployment period was September 2016 to June 2018. Ice drafts were measured on each mooring with an upward-looking acoustic Doppler current profiler. Further details on the instruments associated with the collection of sea ice data, as well as the processing to produce sea ice thickness, are given in Kirillov et al. (2020). Details on the collection of velocity data, as well as temperature and salinity, and their processing are given in Dmitrenko et al. (2020).

Additional sea ice data were provided through the sea ice concentration product OSU-409-a of the EUMETSAT Ocean and Sea Ice Satellite Application Facility (OSI SAF; Lavergne et al., 2010). This product covers 1978–2015 with a daily spatial resolution of 10 km. The breakup date in each grid cell was obtained by identifying the date when sea ice concentration (SIC) dropped below 15% and remained below 15% for three consecutive days. Further details on the processing are given in Kirillov et al. (2020). As discussed in Kirillov et al. (2020), weekly ice charts from the Canadian Ice Service (Tivy et al., 2011) provide a complimentary classification of the ice pack by partial concentration of different stages of development present at each mooring position throughout the winters of 2017 and 2018. Measurements of sea level were taken from hourly tide gauge data collected at the port of Churchill (location shown in Figure 1b).

Beyond the model evaluation in the early studies mentioned in Section 2, we focussed on evaluating the NEMO configuration to be used for the BaySys climate experiments in two ways. As part of the planned BaySys experiments, a historical control experiment was carried out, initially through 2010. This experiment was subsequently extended through 2018 to allow comparison with the BaySys moorings, as well as with additional BaySys sea ice analysis. Given the inclusion of tidal forcing, the model’s tidal fields were compared with the output of the TPXO8 tidal model (Egbert and Erofeeva, 2002) and water level data at the port of Churchill.

Kirillov et al. (2020) examined the sea ice at three HB moorings as well as satellite and CIS data, as detailed in Section 3. Here, we compare the historical control run against the observed ice fields for 2016–2017 and 2017–2018 (Figure 3). For the AN01 mooring, ice formation in the model starts at the same time as in the observations in December 2016. In both, there is a short period of rapid growth to about 0.5 m within a week. Then, there is slow quasi-linear thermodynamic growth through the rest of winter, reaching 1.5 m by April. The model ice is on the thicker end of the observational spectrum and thicker than the empirical line, but not unrealistic. The mooring ice distribution is more variable and the ice is thinner in May than in the model. The model ice thins later than in the observations, but the slope is consistent with the observations. The CIS charts suggest thick first-year ice in July when the model still has about 0.7 m thick ice. The behaviour is similar in 2017, with the model ice at the mooring location forming in November, as with the observations. The sea ice is thinner than in 2016–2017, only reaching about 1.2 m thick by late winter. The model sea ice is thicker than the main distribution but much closer to the empirical line during this year. The model sea ice melts too slowly, with 0.5–1.0 m thick ice remaining in June, before disappearing in July. This result could be partly due to the temperature of the riverine inputs likely being too cold. The NEMO setup used here treats the river water as entering the bay with the same temperature as the grid cell into which it enters, when the terrestrial water inputs likely would be warmer.

The model versus mooring comparison is similar at NE03. In spring 2017, the model sea ice exceeded 2 m in thickness. The discrepancy from the observations is greater during this time period, although some similarly thick ice is in the observed distribution. Decay is faster in the model than at AN01, and the model sea ice disappears by early July, close to the observational timing. The sea ice fields in 2017–2018 also compare well with the model at the upper end of the observed distribution and close to the empirical line. At this location, melting occurs at around the same time as observations for this year, with ice-free conditions in June.

Next, we considered the behaviour at mooring JB02. Here, the modelled sea ice thickness falls in the upper range of observed thickness and is slightly thicker, on average, than the empirical thickness. The model sea ice thickness reaches 1.5 m by May, with a rapid decay in June, albeit slightly slower than observed. The comparison is similar in 2017–2018, when the model sea ice thickness is within the thicker part of the observed distribution. But the model decay in June is much slower than the observations, with July still having 1 m thick sea ice.

We also used the historical control run outputs to evaluate the impact of winter atmospheric forcing on the interannual variability of ice breakups in the HB. To accomplish this evaluation, we applied the same approach that was used in Kirillov et al. (2020) and considered winter as the period when modelled sea ice concentration and reanalysis surface air temperature over central HB (78–95°W, 55–63°N) were greater than 50% and below 0°C, respectively. Then, we plotted the breakup patterns corresponding to the cold years with strong eastward atmospheric transport (type I) and warm years with relatively weak westerly winds (type IV), as discussed in Kirillov et al. (2020; Figure 4). We also calculated the linear correlations between different characteristics of winter atmospheric forcing (freezing degree days, cumulative wind transport and wind vorticity) and modelled breakup dates across Hudson Bay (Figure 5), again similar to Kirillov et al. (2020).

In Kirillov et al. (2020), passive microwave data for sea ice concentration were used to determine the beginning of the winter period over the HB. Here, we used the modelled sea ice concentration that led to slightly different lengths and the cumulative effect of atmospheric forcing during winter from 1980 to 2018 (see the differences in our Figure 4a and figure 6a in Kirillov et al., 2020).

The model reproduces the general structure of the mean breakup patterns with either type of atmospheric forcing, with the latest breakup in southern HB and the earliest breakup in the northwest corner. In general, the model breakup dates are late by about a month compared to the observations. Yet, the differences between the sea ice retreat dates for type I and IV forcings are similar to the observations. Thus, the largest difference of up to 30–40 days are found in the southeastern part of the bay, while the smallest differences are localized in the northwest, which coincides with the pattern found from satellite measurements. The model reproduces the impact of changes in atmospheric forcing during winter on the summer sea ice retreat adequately: the linear correlations between different atmospheric impact factors and the breakup dates compare well with those found from the observational data (our Figure 5 versus figure 7 in Kirilov et al., 2020).

These sea ice production patterns were also observed in the evaluation of the biogeochemical model, BLINGv0 + DIC, in Deschepper et al. (2023) for 2018. The delay in the melt and production of sea ice impacts the timing of the spring and autumn blooms of chl-a in HBC. However, the spatial patterns observed in the models were similar to those seen by MODIS-Aqua satellite chl-a data in 2018. Further evaluation of the models showed that the variability of spatial and temporal patterns was influenced by the mixed layer depth dynamics (18% to 14% of variability), river runoff (13% to 6% of variability) and sea ice production (7% of variability) for the MIROC5 and MRI forced simulations, respectively (Deschepper et al., 2023). These findings illustrate the usefulness of physical and biogeochemical oceanographic models in understanding the impacts that physical drivers may have on the base of the food web.

Here, we compare the model temperature, salinity and current velocity fields with the same three BaySys moorings (Figure 1b) at the various depths where the sensors were mounted. To avoid high-frequency fluctuations, we compare the monthly means for each data set. Because the model outputs 5-day means, the variability within a month will be smaller compared to the higher frequency observations.

At AN01, temperature and salinity were measured at 36 m, 55 m and 103 m. At 36 m and 55 m, both observational and model temperatures remain close to the freezing point in winter (January through June, Figure 6a and b), when sea ice is present at the surface (Figure 3). The model shows greater warming at these two levels during the ice-free season relative to the observations, with differences reaching up to 2°C in October of each year. We suspect that shortwave radiation penetrates too deeply in the model in summer, which then leads to the warm bias in Autumn. At 103 m, the temperature range is more constrained, although the amplitude of the model seasonal cycle is larger than in the observations. For the salinity, at both 36 m and 55 m, the model and observations track very well in 2018. The discrepancies are larger in 2017, and more noticeable at 36 m than 55 m. The model salinity at 103 m behaves similarly to 55 m, with a freshening event in winter 2017–2018 that is not seen at the mooring.

At NE03, the model and observational comparison is very good (Figure 7). Here, both time series show the same near-freezing temperatures at both 28 m and 43 m throughout the ice season. As well, both model and observations capture the summer warming at both depths, with both reaching a maximum of 4–5°C in October 2016. That said, the model salinity is fresher than the mooring measurements at both depths in 2017, but the difference becomes significantly smaller in 2018. Abnormally high discharge conditions existed from 2016 to 2018, which may create extended freshwater plumes at different locations throughout the HBC. Also, there is much more noise in the ice data in 2017–2018, which may impact the salinity. Additionally, 2017 was a record-high flow year in the Nelson and Churchill rivers.

At JB02, the model and observational temperatures compare well during 2016–2017 at all depths, including for the warm period in the autumn of 2016 (Figure 8). The model warms faster than the observations in late summer 2017 when the observational record unfortunately stops. The model salinity is slightly larger than the mooring salinities in winter 2017 at 35 m and 54 m, although the two records overlap nicely at 74 m. During the last part of the observational record, in the spring and summer of 2017, the model salinities fall within the range of the observed salinities at all depths.

Current rose diagrams are used to compare velocity magnitude and directions at each mooring and instrumented depth (Figure 9). At AN01, both model and observations show, at 10 m and 20 m, that the mean flow is to the southwest, with similar magnitudes, although the model does not capture the largest values in the observational record. Periods when the flow is to the north and east in the observational record are also captured at these depths, even if the likelihood of these events is underestimated. The model flow remains more barotropic than the observations (a known issue with ocean models; Penduff et al., 2005). Thus, aspects of the flow variability observed at 40 m are not captured in the model.

Similarly, at NE03 (Figure 10), the main eastward component of the flow from the mooring record is captured in the model. However, the variability at depth, especially near the bottom at 40 m, is not captured. At JB02 (Figure 11), the model captures the northward flow out of James Bay at all depths in terms of both direction and magnitude. Some of the variability in flow direction at 20 m is not reproduced, and the model currents at 40 m are slower than the observations. Overall, the current rose diagrams demonstrate a generally good correspondence between the predominant current directions and speeds in the modelled and observed data sets.

In terms of examining the representation of tides in the ANHA configuration, we compared the model solution with that from TPXO8 (Egbert and Erofeeva, 2002). We considered the four constituents with the largest amplitudes in the bay: M2, S1, N2 and K1. For each constituent, we chose 500 random points to compare amplitude and phase (Figure 12). For constituents like K1, the model and TPXO8 amplitudes are nearly equivalent, with a slight overestimation of the amplitude in the ANHA configuration. For the other constituents, the ANHA configuration consistently underestimates the amplitude by a factor of 2–4. This underestimation is especially true in Hudson Strait and Ungava Bay, where the M2 tides are known to be very large. This outcome is not surprising considering the lower resolution in the ANHA4 configuration used here (Kleptsova and Pietrzak, 2018). Additionally, the bathymetry of the HBC is known to have many uncertainties, impacting modelled amplitudes. However, the phases agree very well for both products, for all constituents.

Additionally, we compared model estimates of sea level with observed tide gauge anomalies at the Port of Churchill (Figure 13). Both 5-day averages for the period between 2004 and 2015 and the mean monthly averages over the same time period show a good correspondence between the model and observations. The correlation coefficient between these two records is 0.949, significant at the 99% level.

The objective of the modelling team in BaySys was to support the other BaySys teams in investigating the relative impacts of climate change and regulation on freshwater-marine coupling within the HBC (Foxe Basin, Hudson Bay, James Bay, Hudson Strait and Ungava Bay). In support of the Team 1 (Physical Oceanography) hypotheses, a sea ice and oceanographic model was used to further study the effects of freshwater loading and ice cover on the circulation of HB. This modelling perspective is based on the NEMO ocean general circulation model coupled to the LIM2 sea ice model. Central also to Team 6 (Modelling) goals was the development of an integrated observational-modelling framework that will provide insight into and improved representation of the physical, biological and biogeochemical processes in the Hudson Bay system. The modelling will provide a framework and tool with which to simulate projected changes in marine state and dynamic variables, while also enabling an integration of observations and numerical analyses.

Team 6 thus focused on the application of a modelling framework for the BaySys project that will provide insight into the relative effects of climate change and hydroelectric regulation on physical and biogeochemical conditions in the HB system. Thus, an existing NEMO modelling configuration, ANHA4, was selected for the BaySys project. The version of the ANHA NEMO configuration that existed at the start of the project was used in a number of initial studies, discussed in greater detail in Section 2, that helped provide the framework for the developments needed to carry out the planned long BaySys climate change integrations. Such developments include switching to the newer v3.6 of NEMO that allowed for the inclusion of explicit representation of the tides and also enhanced the speed of the simulations through the inclusion of land masking (Madec et al., 2008). Significant effort was spent in incorporating runoff from the new and improved HYPE hydrological models (Stadnyk et al., 2021), as well as incorporating forcing from the bias-corrected versions of multiple CMIP5 models. Improvement also occurred in terms of the biogeochemical modules, with the inclusion of a carbon module as part of BLING, and the ongoing work developing and coupling BIGCIIM, with its sympagic as well as pelagic components. In the end, the NEMO ANHA4 configuration, in conjunction with the BLINGv0 + DIC model and seven passive tracers, took only approximately 18 hours of CPU time when run on 256 processors on the Digital Research Alliance, Canada, machine graham, to run 2 years of simulation. This efficiency made it possible to carry out the 10 near-century-long integrations (5 with naturalized runoff and 5 with regulated runoff) central to the BaySys goal of studying the relative impacts of climate change and regulation. Yet, in terms of real time, admittedly dependent on external factors like throughput on the Digital Research Alliance systems, which is driven by allocation size, about 3–4 months of running time was required per experiment. This requirement clarifies why the NEMO model ensembles have only five members for each runoff case, while the hydrological model used 19 (Stadnyk et al., 2021). Braun et al. (2021) demonstrated that the 5 NEMO climate simulations represented the ensemble of 19 across a multitude of climatic indicators reasonably well. A larger NEMO model ensemble would have increased the reliability of the predicted outcomes, but was not practical computationally.

As part of any modelling study, the results need to be evaluated against observations and other models to understand the model’s ability to simulate the real world properly and understand the limitations and weaknesses of the model. In general, because no model is good at representing all processes and all scales, such evaluation has to be carried out as part of the detailed analysis within individual studies. Thus, aspects of this work within BaySys are being carried out and reported in each of the modelling studies using NEMO output, some of which are highlighted below.

Also suggested is a need for greater interactions between modellers and observationalists early in the planning processes to ensure that the observations are appropriate for evaluating the model, especially the biogeochemical elements of it. Although not discussed in this manuscript, inter-model comparison was also carried out during the BaySys project. However, it soon became apparent that trying to compare the output from a short-range prediction system (e.g., Regional Ice Ocean Prediction System, RIOPS, of the Meteorological Service of Canada) with extended free-running prognostic simulations for climate applications is very challenging due to the very different timescales involved. Again, the key message is that such comparisons need to be built into projects from the start, and both time and people need to be included to allow such comparisons to have significant value to a project and lead to model improvement. Sufficient time has to be provided to ensure that knowledge gained from model-observation and model-model comparisons can be brought back into a modelling system, to allow schemes and parameters to be fully tested before long climate scenarios are run.

That said, some general overall big-picture evaluation was carried out in Section 4 to compare the historical control run with the BaySys moorings, as well as a subset of other observations available within the bay. The historical control run was used for this evaluation given that it is based on realistic forcing from a reanalysis product (ERA-Interim), unlike the climate experiments, which are forced from given climate simulations with their own internal variability. As the general evaluation in this article shows, in many ways, the model agrees well with the observations, reproducing well the observed seasonality of temperature, salinity, and sea ice, as well as the tidal cycle, especially when considering the inherent limitations of comparing point source measurements with model fields that are averaged over a grid cell. Given that general agreement, there are discrepancies between model and observations, some of which suggest too much shortwave heat penetration in summer, which will need to be considered carefully when discussing model results. Overall, the authors are satisfied with the model configuration that has been run for BaySys, representing that it does well finding the main features of the circulation and hydrography of the HBC. Thus, the BaySys model experiments provide an appropriate basis for beginning to understand freshwater-marine coupling and the relative role of climate change and regulation in the HBC.

Results using output from these experiments will appear in many studies, both in this special feature and beyond. Dmitrenko et al. (2020) used model output to help scale up mooring data to the bay scale to show that atmospheric vorticity sets the basin-scale circulation within Hudson Bay. Central to BaySys is the question of the relative role of climate change versus river regulation. This question led to two core papers putting the present day in the context of the longer historical record (Lukovich et al., 2021b; 2021c) before using the long integrations to show that climate change impacts are evident in terms of sea surface temperature increases and sea ice decreases (Lukovich et al., 2021a). The latter (Lukovich et al., 2021a) also showed that regulation suppresses climate impacts in winter and reinforces them in summer. Additionally, an analysis of the long-term trends in bay-wide sea ice, circulation and hydrography will be provided in future articles. NEMO results are also being used to help understand the impact of improved representation of river runoff through the A-HYPE hydrological model in Stadnyk et al. (2021). For the biogeochemical analysis, NEMO output is being used to drive a box model as well as to study how the biogeochemical processes within the bay are represented by BLING.

As the large numbers of authors indicate, getting a suitable NEMO model configuration developed and running for the BaySys project was a major endeavour. A huge amount of development, testing and evaluation occurred at intermediate stages that are not highlighted in any publication, but was needed to produce the present product. In the end, we believe that such an effort was well worth it, producing a tool that allows for the detailed study of the HBC, now and in the future, for many years to come. As part of this process, more experiments and more years of simulation of the oceanographic conditions within the HBC have been carried out compared to all the previous modelling studies discussed in Section 2 combined.

Continuing work and development on the model and the modelling system is warranted. Gaining access to improved bathymetry is likely to improve the representation of the circulation. Examining how shortwave radiation penetrates the upper water column will likely improve the representation of the upper ocean structure. Switching to a more advanced multi-category sea ice model would potentially improve issues, especially regarding breakup and freeze-up dates. Resolution could be enhanced (at the cost of significantly increased computational cost) to look at current structure and coastal processes in more detail. Potentially, the use of nesting software may allow for detailed studies of given regions, such as estuaries. In the end, although comprehensive, the modelling in BaySys is just a start for simulating the circulation and hydrography, and their evolution, in the HBC.

Mapped fields of total alkalinity and dissolved inorganic carbon are from GLODAPv2. These data were prepared by the Carbon Dioxide Information Analysis Center for the Climate Change Research Division Office of Biological and Environmental Research, U.S. Department of Energy. These data are freely available at https://www.ncei.noaa.gov/data/oceans/ncei/ocads/data/0162565/. Mapped fields of objectively analysed macronutrient (nitrate and phosphate) and oxygen climatology are freely available through the National Centers for Environmental Information (NCEI https://www.nodc.noaa.gov/OC5/woa13/woa13data.html), which is funded by the National Oceanic and Atmospheric Administration (NOAA). The CESM1 project is supported primarily by the National Science Foundation. The material is based upon work supported by the National Center for Atmospheric Research, which is a major facility sponsored by the National Science Foundation under Cooperative Agreement No. 1852977. We thank all the scientists, software engineers, and administrators who contributed to the development of the CESM1. Data can be accessed at https://www.earthsystemgrid.org/search.html?Project=CMIP5&Ensemble=r1i1p1&Model=CESM1-BGC&Product=output1. All simulation output used in this research is available on Compute Canada (www.computecanada.ca) and can be requested through the Canadian NEMO Forum (https://canadian-nemo-ocean-modelling-forum-commuity-of-practice.readthedocs.io/en/latest/index.html).

This work is a contribution to the Natural Sciences and Engineering Council of Canada (NSERC) Collaborative Research and Development project titled BaySys. In addition, this research contributes to the ArcticNet and MEOPAR Networks of Centers of Excellence and the Arctic Science Partnership.

Funding for this research was graciously provided by Manitoba Hydro, NSERC, Amundsen Science, and the Canada Research Chairs program.

The authors have no competing interests, as defined by Elementa, that might be perceived to influence the research presented in this manuscript. Jean-Éric Tremblay is an Associate Editor at Elementa. He was not involved in the review process of this article.

Contributed to conception and design: PGM, DB, MB, JL, TP, TAS, KS, JET.

Contributed to acquisition of data: PGM, LB, LCG, ID, YG-Q, LCG, NG, XH, JMM, CP, NR, TAS, RT, AT, YX.

Contributed to analysis and interpretation of data: PGM, LB, LCG, ID, YG-Q, LCG, XH, SAK, SJ, JL, JMM, CP, NR, TAS, RT, AT, YX.

Drafted and/or revised the article: PGM with assistance of all other articles.

Approved the submitted version for publication: All authors.

How to cite this article: Myers, PG, Barber, D, Braun, M, Buchart, L, Castro de la Guardia, L, Deschepper, I, Dupont, F, Ehn, J, Garcia-Quintana, Y, Gillard, LC, Grivault, N, Hu, X, Kirillov, SA, Jafarikhasragh, S, Lukovich, J, Maps, F, Marson, JM, Papakyriakou, T, Pennelly, C, Ridenour, N, Stadnyk, TA, Sydor, K, Tao, R, Tefs, A, Tremblay, J-É, Xu, Y. 2024. An overview of the NEMO modelling for the BaySys project. Elementa: Science of the Anthropocene 12(1). DOI: https://doi.org/10.1525/elementa.2022.00111

Domain Editor-in-Chief: Jody W. Deming, University of Washington, Seattle, WA, USA

Knowledge Domain: Ocean Science

Part of an Elementa Special Feature: The Hudson Bay System Study (BaySys)