A regional physical–biogeochemical ocean model for marine resource applications in the Northeast Pacific (MOM6-COBALT-NEP10k v1.0)

Regional ocean models enable the generation of computationally affordable and regionally tailored ensembles of near-term forecasts and long-term projections of sufficient resolution to serve marine resource management. Climate change, however, has created marine resource challenges, such as shifting stock distributions, that cut across domestic and international management boundaries and have pushed regional modeling efforts toward “coastwide” approaches. Here, we present and evaluate a multidecadal hindcast with a Northeast Pacific regional implementation of the Modular Ocean Model, version 6, with sea ice and biogeochemistry that extends from the Chukchi Sea to the Baja California Peninsula at 10 km horizontal resolution (MOM6-COBALT-NEP10k, or NEP10k). This domain includes an Arctic-adjacent system with a broad, shallow shelf seasonally covered by sea ice (the eastern Bering Sea), a sub-Arctic system with upwelling in the Alaska Gyre and predominant downwelling winds and large freshwater forcing along the coast (the Gulf of Alaska), and a temperate, eastern boundary upwelling ecosystem (the California Current Ecosystem). The coastwide model was able to recreate seasonal and cross-ecosystem contrasts in numerous ecosystem-critical properties including temperature, salinity, inorganic nutrients, oxygen, carbonate saturation states, and chlorophyll. Spatial consistency between modeled quantities and observations generally extended to plankton ecosystems, though small to moderate biases were also apparent. Fidelity with observed zooplankton biomass, for example, was limited to first-order seasonal and cross-system contrasts. Temporally, simulated monthly surface and bottom temperature anomalies in coastal regions (<500 m deep) closely matched estimates from data-assimilative ocean reanalyses. Performance, however, was reduced in some nearshore regions coarsely resolved by the model's 10 km resolution grid and for point measurements. The time series of satellite-based chlorophyll anomaly estimates proved more difficult to match than temperature. System-specific ecosystem indicators were also assessed. In the eastern Bering Sea, NEP10k robustly matched observed variations, including recent large declines, in the area of the summer bottom water “cold pool” (<2 °C), which exerts a profound influence on eastern Bering Sea fisheries. In the Gulf of Alaska, the simulation captured patterns of sea surface height variability and variations in thermal, oxygen, and acidification risk associated with local modes of interannual to decadal climate variability. In the California Current Ecosystem, the simulation robustly captured variations in upwelling indices and coastal water masses, though discrepancies in the latter were evident in the Southern California Bight. Enhanced model resolution may reduce such discrepancies, but any benefits must be carefully weighed against computational costs given the intended use of this system for ensemble predictions and projections. Meanwhile, the demonstrated NEP10k skill level herein, particularly in recreating cross-ecosystem contrasts and the time variation of ecosystem indicators over multiple decades, suggests considerable immediate utility for coastwide retrospective and predictive applications.

1 Introduction

The western coasts of the continental United States, Canada, and Mexico form the eastern bounds of the North Pacific Gyre, which substantially impacts North American climate and supports a diverse assemblage of ecosystems, species, and resources. These ecosystems include valuable fisheries that represented roughly 42 % of the USD 4.6 billion in commercial U.S. domestic landings in 2020 (National Marine Fisheries Service, 2022). Management of these interconnected, multiscale marine resources presents a challenge, particularly with the growing need to account for changing climate and ocean conditions. Ocean warming, acidification, and deoxygenation stand to fundamentally alter coastal ecosystems (Gruber, 2011), potentially driving fluctuations in living marine resource abundance due to habitat range shifts (e.g., Pinsky et al., 2013; Christian and Holmes, 2016; Smith et al., 2021; Chasco et al., 2022; Thompson et al., 2023), recruitment and fish size changes (e.g., Holsman et al., 2019; Litzow et al., 2022), and heightened competition and predation from invasive species (Grosholz et al., 2000; Zeidberg and Robison, 2007; Compton et al., 2010). Additionally, extreme events such as marine heat waves (e.g., Rogers-Bennett and Catton, 2019; McPherson et al., 2021) and harmful algal blooms (e.g., Anderson et al., 2015) can degrade foundational habitats and compromise water quality.

Numerical ocean models facilitate both the understanding of difficult-to-observe ocean and ecosystem dynamics and the forecasting and projection of near- to long-term ocean conditions. Previous regional modeling efforts in the Northeast Pacific Ocean have contributed considerably to our understanding of the Bering Sea (Danielson et al., 2011; Hermann et al., 2013, 2016; Cheng et al., 2015; Pilcher et al., 2019; Kearney et al., 2020), Gulf of Alaska (Hermann et al., 2009; Hinckley et al., 2009; Cheng et al., 2012; Coyle et al., 2012, 2019; Danielson et al., 2020; Hauri et al., 2020, 2024), and California Current System (Marchesiello et al., 2001; Di Lorenzo et al., 2005; Gruber et al., 2006; Veneziani et al., 2009; Neveu et al., 2016; Van Oostende et al., 2018; Dussin et al., 2019; Deutsch et al., 2021; Renault et al., 2021) and the broader NEP10k domain (Desmet et al., 2022, 2023). Predictions and projections from these regionally tailored ocean models have also been enlisted to understand and anticipate living marine resource responses to climate variability and change (e.g., Gruber et al., 2012; Hermann et al., 2016; Holsman et al., 2020; Siedlecki et al., 2016; Howard et al., 2020; Pozo Buil et al., 2021; Pilcher et al., 2022; Jacox et al., 2023). In a growing number of cases, applications have been extended to management (e.g., Anderson et al., 2016; Punt et al., 2021; Brodie et al., 2023; Smith et al., 2023; Hollowed et al., 2024). Such applications have been hampered, however, by the use of relatively small domains and limited ensembles to characterize uncertainties. Climate change impacts and species responses traverse the bounds of those domains, thus motivating an integrated “coastwide” modeling framework with rigorously defined uncertainties.

A key challenge is thus configuring a coastwide modeling framework with sufficient resolution and complexity to adequately represent fisheries-critical ocean features across the full domain while also maintaining low computational cost conducive to generating ensembles (Drenkard et al., 2021). This challenge is made more acute by the diversity of Northeast Pacific ecosystems and the mechanisms by which climate shapes them. The Bering Sea, for example, features one of the world's broadest shallow continental shelf environments, which supports benthic and demersal fisheries that are amongst the most productive in the world (National Research Council, 1996). These fisheries, however, have proven to be highly sensitive to temperature and food fluctuations in these shallow habitats (Hunt et al., 2002, 2011). Recent warming and reduced food supply in the eastern Bering Sea (EBS), for example, was linked to the collapse of the snow crab fishery (Szuwalski et al., 2023). Productivity as well as benthic and pelagic habitat fluctuations on the eastern Bering shelf are further linked to coupled ocean and sea ice dynamics (Mueter and Litzow, 2008; Brown and Arrigo, 2013; Hunt et al., 2022), presenting an additional challenge for ocean modeling systems intended for fisheries applications in this region.

In the Gulf of Alaska (GOA), downwelling winds and abundant freshwater input prevail and contribute to a strong cyclonic circulation of the Alaska Gyre (Stabeno et al., 2004). Despite the predominance of downwelling winds, the confluence of the high nitrate waters of the basin with the high iron waters of the shelf (assisted by shelf-break eddies), as well as upwelling of nitrate by wind stress curl, promote high production in the coastal GOA (Stabeno et al., 2004; Hermann et al., 2009; Coyle et al., 2019). While correlation with the El Niño–Southern Oscillation (ENSO) can be found (e.g., Bailey et al., 1995; Whitney and Welch, 2002; Amaya et al., 2023b), lower-frequency modes of decadal climate variability tend to predominate (e.g., Di Lorenzo et al., 2008) and are associated with marked decadal-scale ecosystem regime shifts (Anderson and Piatt, 1999; Hare and Mantua, 2000) and modulations in fisheries and ecosystem risks (Hauri et al., 2021, 2024). Cold water temperatures and the proximity of North Pacific basin waters, which are exceptionally rich in dissolved inorganic carbon (DIC), make the Gulf of Alaska particularly susceptible to ocean acidification (Fabry et al., 2009; Byrne et al., 2010; Mathis et al., 2015). Periodic on-shelf intrusions of DIC-rich deep Pacific water can suppress the aragonite and calcite saturation states and stress commercially important crab and shell fisheries (Ladd et al., 2005). Increased freshwater input due to deglaciation, which is naturally low in alkalinity, may also exacerbate coastal acidification trends (Reisdorph and Mathis, 2014; Evans et al., 2014). In offshore waters, the iron supply strongly modulates ocean productivity, though the impacts of such variations on fisheries remains speculative (Lippiatt et al., 2010; McKinnell, 2013; Kearney et al., 2015).

The California Current is one of the four major eastern boundary upwelling systems in the global ocean (Hill et al., 1998). Marine resource fluctuations are inextricably linked to variations in the timing, strength, and source waters of this seasonal upwelling (e.g., Bograd et al., 2009). Physical, biogeochemical, and marine resource dynamics of the California Current correspond strongly with ENSO (Ohman et al., 2017; Turi et al., 2018; Cordero-Quirós et al., 2022) through diverse atmospheric and oceanic teleconnection pathways (Alexander et al., 2002; Jacox et al., 2015; Frischknecht et al., 2015). While a narrow shelf and modest riverine inputs over much of the coast give the California Current an oceanic character, the system nonetheless supports significant benthic and demersal fisheries, which are periodically subject to heightened hypoxia and acidification risks common in upwelling systems (Bograd et al., 2008; Hauri et al., 2009; Wolfe et al., 2023). These risks can be further amplified by processes resulting from changing land use such as increased nutrient input, pollution, and coastal engineering (e.g., Halpern et al., 2009; Hughes et al., 2015). The considerable productivity generated by coastal upwelling also supports climate-sensitive forage fish, highly migratory species, and top predators that are ecologically, economically, and culturally important. Projections suggest that upwelling strength, seasonality, and source water properties may shift with climate change (Rykaczewski and Dunne, 2010; Rykaczewski et al., 2015; Sydeman et al., 2011; Pozo Buil et al., 2021) and significantly alter ecosystem productivity and fisheries (McClatchie et al., 2010; Bograd et al., 2023; Jacox et al., 2024).

Here, we present a regional implementation of the modular ocean model (MOM6) with coupled sea ice and biogeochemistry spanning the Northeast Pacific and assess the degree to which this system can capture fisheries-critical mean patterns and fluctuations across the diverse ecosystems of the Northeast Pacific. We evaluate the model's capacity to represent both large-scale contrasts in ecologically important variables across ecosystems and variations in fisheries-oriented diagnostics within each ecosystem. We also assess computational costs to ensure the feasibility of ensemble predictions. We conclude with an assessment of the model's current utility for fisheries applications and a discussion of priority developments for addressing model biases in order to maximize future utility in informing fisheries and ecosystem decisions.

2 Methods

2.1 Physical model configuration

The NEP10k model domain (Fig. 1) is designed to cover the western coast of the continental United States and contiguous regions. It extends from 10.8–80.7° N and 156.6° E–105.0° W, measuring 3320 ± 126 km by 7764 ± 58 km (mean ± standard deviation) in the offshore and along-shore dimensions, respectively. The model is integrated on an orthogonal curvilinear grid that consists of 342×816 tracer cells with a horizontal resolution averaging 9.7 ± 0.5 km and a minimum bathymetric depth of 10 m. The domain has four open boundaries, the longest of which arcs through the Pacific Ocean and is referenced as the “western” boundary. In the vertical direction, the model uses 75 z* coordinates, which are approximately consistent with the depth-from-mean-sea-level but are stretched by variations in sea surface height across all water column layer thicknesses rather than isolating that variability in the surface layer (Adcroft and Campin, 2004). We prescribe a layer thickness of 2 m from the surface to a depth of 8 m, between 2.01 to 2.34 m in thickness between 8 and ∼31 m in depth, then with spacing gradually increasing to 250 m in the deepest portions of the model domain. Bathymetry for the NEP10k domain was derived from the 2020 General Bathymetric Chart of the Oceans (GEBCO Bathymetric Compilation Group, 2020) and is not vertically rounded or truncated. MOM6 does not need the topography to conform to the vertical level thicknesses but instead can let the bottommost non-vanished layer vary in thickness to match the topography and then collapse the layer to zero thickness below where the topography intersects the model layer. Simulations used a baroclinic time step of 400 s and a variable barotropic time step set to maintain stability (Hallberg, 1997; Hallberg and Adcroft, 2009). A longer, 1200 s time step was used for thermodynamic and biogeochemical tracer calculations, as thermodynamic processes tend to evolve more slowly than dynamic ones. Past studies have used a longer time step for these processes without compromising their representation while reducing the overall computation time (e.g., Ross et al., 2023). The success of this strategy for the NEP10k domain will be assessed herein.

Figure 1 NEP10k domain and bathymetry. NEP10k domain and bathymetry with a log-normal color scale to emphasize priority coastal regions. White coloration indicates non-ocean (i.e., land-masked) grid cells that are not computed in model integrations, which include the Sea of Okhotsk. The agglomerate land mask is outlined in black. Red lines indicate the areas that are spatially averaged for regional shelf temperature and chlorophyll time series. These regions, from north to south, are the Bering Sea (BS), Gulf of Alaska (GOA), British Columbia (BC), Northern California Current System (NCCS), Central California Current System (CCCS), and Southern California Current System (SCCS). The southern arc of the Bering Sea polygon traces the Aleutian island chain; the southernmost land bounds of the Southern California Current System and Gulf of Alaska polygons, as well as both the northernmost and southernmost land bounds of the British Columbia polygons, roughly correspond with international geopolitical boundaries. The dark green contour delineates the 500 m isobath, which we use to isolate shelf grid cells (i.e., where depth ≤500 m).

The core components of the physical ocean model, Modular Ocean Model 6 (MOM6), are described in Adcroft et al. (2019). A full account of the parameterization choices implemented for the simulations presented in this study can be found in the Supplement in Drenkard et al. (2024a) (MOM_parameter_doc.all). Here, we elaborate on a few choices (Table A1), highlighting consistencies and contrasts with the recently published Northwest Atlantic configuration documented in Ross et al. (2023). As in Ross et al. (2023), ocean boundary layer mixing, specifically vertical turbulent mixing coefficients in the surface layer, are parameterized using the energetic planetary boundary layer (ePBL) scheme developed by Reichl and Hallberg (2018). However, unlike Ross et al. (2023), we switched to the submesoscale mixing and restratification scheme of Bodner et al. (2023) from that of Fox-Kemper et al. (2011). The Bodner parameterization has the advantage of dynamically calculating the submesoscale front length (i.e., the length scale perpendicular to the front), which can vary significantly seasonally and latitudinally across the ecosystems represented in NEP10k (Bodner et al., 2023). In the ocean interior below the surface boundary layer, mixing primarily depends on the shear-driven turbulence mixing scheme of Jackson et al. (2008). The standard Jackson formulation, however, was found to overmix some shelf regions subject to strong tidal motions. This overmixing was ameliorated by including a scaling factor for the turbulent decay length scale. Bottom drag and horizontal viscosities were parameterized as in Ross et al. (2023). Unlike Ross et al. (2023), the background kinematic viscosity parameter, KV, was set to 0.0 m² s⁻¹; this parameter is intended to supplement the existing dynamic viscosity (based on the diapycnal diffusivity, KD) and was determined to be unnecessary for this application. Sea ice is modeled with the Sea Ice Simulator, version 2 (SIS2, Adcroft et al., 2019). This sea ice model uses five sea ice thickness categories and no explicit ridging scheme. The sea ice rheology is an elastic-viscous-plastic scheme (Hibler, 1979) and a directionally split piecewise constant advection scheme for thickness. The delta-Eddington radiation scheme is used, and the internal thermodynamics are enthalpy-conserving (Briegleb and Light, 2007).

2.2 Physical model forcing

The ocean hindcast simulation was run from 1993 through 2019 on NOAA's GAEA supercomputer, which is housed and managed in partnership with the Department of Energy through the National Climate-Computing Research Center. Hourly atmospheric forcing for NEP10k was prescribed from the European Centre for Medium-Range Weather Forecasts Reanalysis 5 (ERA5; Hersbach et al., 2020). The bulk formulae of Large and Yeager (2004) were used to calculate latent and sensible heating after adjusting to the 2 m ERA5 reference height. Light attenuation and associated heating within the water column was calculated from Manizza et al. (2005) using dynamically varying chlorophyll from the biogeochemical model (Sect. 2.3).

Daily freshwater runoff was prescribed using output from the Global Flood Awareness System, version 4.0 (GloFAS; Harrigan et al., 2020; Grimaldi et al., 2022) – a hydrological inundation model that is also forced by ERA5. Freshwater discharge at ocean-adjacent “pit cells” in GloFAS was remapped to the nearest MOM6 coastal ocean grid cells. Pit cells are GloFAS grid cells where the local drain direction indicates that only inward water flow occurs and is therefore a point of accumulation (e.g., lakes) or a point of egress to the ocean via either ocean adjacency or connectivity through other pit cells (e.g., wetlands). For the Gulf of Alaska, we substituted freshwater discharge from Beamer et al. (2016; Hill, 2023), a model dedicated to the representation of freshwater discharge and glacier mass balance in Alaska, with calibration against observed watersheds.

Open lateral boundary and initial conditions for temperature, salinity, sea surface height, and momentum were prescribed as daily means from the $\mathrm{1}/\mathrm{12}$° Global Ocean Physics Reanalysis (GLORYS12; Lellouche et al., 2021). Tidal forcing was prescribed at the boundaries using the amplitude and phase from the global tidal elevation and transport atlas, version 9 (TPXO; Egbert and Erofeeva, 2002). Tides were implemented as in Ross et al. (2023) with four semidiurnal constituents (M2, S2, N2, K2), four diurnal constituents (K1, O1, P1, Q1), and two long-period constituents (Mm and Mf). Initial and boundary conditions were regridded to the NEP10k domain using the xESMF Python software package (Zhuang et al., 2023). Boundary conditions were imposed as in Ross et al. (2023), with barotropic flows handled with a Flather (1976) boundary condition, while baroclinic flows are handled with an Orlanski (1976) radiation condition. Lateral boundary forcing also applies nudging and tracer reservoirs (the latter retains a memory of water properties exchanged with the modeling domain rather than instantaneous forcing; see Ross et al., 2023, for more details). As in Ross et al. (2023), the lateral ocean boundary radiation and nudging schemes utilize 3 d inflow, 360 d outflow timescales, and both inward and outward tracer reservoir length scales were 9000 m (Table A1). No nudging was included in the interior of the domain.

2.3 Biogeochemical model configuration

Biogeochemistry was simulated using version 3.0 of the Carbon, Ocean Biogeochemistry, and Lower Trophics (COBALTv3.0) model (Stock et al., 2025; Ross et al., 2023). COBALTv3.0 includes 40 prognostic state variables to capture plankton food web dynamics and the cycling of carbon, nitrogen, phosphorus, iron, silica, calcium carbonate, and lithogenic material in ocean and coastal environments. COBALTv3.0 builds on prior COBALT formulations (Stock et al., 2014, 2020) by adding a third phytoplankton size class following Van Oostende et al. (2018). The resulting small, medium, and large sizes correspond to the canonical pico-, nano-, and microplankton size classes defined by Sieburth et al. (1978) and enable COBALT to better resolve the range of phytoplankton communities from oligotrophic gyres to intensely productive upwelling systems. These join diazotrophs to give a total of four phytoplankton functional types to go along with a plankton food web including three zooplankton functional types and free living bacteria (Stock et al., 2014, 2020). Additional flexibility in zooplankton feeding, direct phytoplankton sinking, and improved photoadaptation and photoacclimation dynamics were also added (Stock et al., 2025), and the formulation enlists an adaptation of the dynamic N:P ratio scheme proposed by Galbraith and Martiny (2015) and initially presented in Ross et al. (2023).

Initial and boundary conditions for biogeochemistry were drawn from the same sources as Ross et al. (2023). The 2018 World Ocean Atlas (WOA18) was used for macronutrients (NO₃, PO₄, SiO₄) and oxygen (O₂), with seasonal averages above 800 m and annual climatologies below (Boyer et al., 2019; García et al., 2019a, b). The Empirical Seawater Property Estimation Routines Locally Interpolated Regressions (ESPER_LIR) presented by Carter et al. (2021) were used to provide initial and time-varying (i.e., seasonal, interannual to decadal variability, and multidecadal trends) boundary conditions for dissolved inorganic carbon and alkalinity. The input values used for this calculation were the location, temperature, salinity, and date. Boundary conditions for other tracers, which generally come into more rapid equilibrium with interior conditions, were drawn from an earlier global ocean hindcast (Stock et al., 2014).

River carbon, alkalinity, nutrients (N, P, and Si) and oxygen inputs were derived by combining the River Chemistry for U.S. Coast (RC4USCoast) database (Gomez et al., 2023) for U.S. waters in the continental United States, the Global River Chemistry database (GLORICH, Hartmann et al., 2019, 2014) for subarctic/Canadian waters, and the Arctic Great Rivers Observatory (Holmes et al., 2012; ArcticGro, 2024). To force COBALT, riverine nutrient inputs are needed for dissolved inorganic and organic nitrogen and phosphorus, particulate nitrogen, phosphorus, and iron. Direct information on dissolved and particulate organic nutrient inputs was not available in all cases. In cases where one or both of these values were missing, the ratio of dissolved and/or particulate organic inputs to dissolved inorganic nitrogen was estimated from the Global Nutrient Export from WaterSheds (Global NEWS) model (Mayorga et al., 2010). This NEWS-derived ratio was then multiplied by the observed inorganic nitrogen to estimate dissolved and particulate organic fluxes in a manner that preserved their relative importance but avoided regional biases in global nutrient-load models such as Global NEWS. Dissolved organic nitrogen and phosphorus was partitioned into 40 % labile, 30 % semi-labile, and 30 % semi-refractory components in COBALT to be consistent with mean tendencies reported by Wiegner et al. (2006). Particulate phosphate is often the largest phosphorus source in rivers, but much of it is buried in nearshore waters before reaching the ocean. Following Froelich (1988), we assumed that 30 % of the particulate phosphorus was mobilized in estuarine sediments to phosphate, with the rest buried. Iron concentrations for all rivers were set to 70 nM (de Baar and de Jong, 2001). As in Ross et al. (2023), atmospheric CO₂ was set using the monthly historical time series of Meinshausen et al. (2017), updated after 2014 using SSP2-4.5 scenario values (Meinshausen et al., 2020), and nutrient, dust, and iron deposition were based on a 1993–2014 climatology from GFDL's ESM4.1 model (Dunne et al., 2020; Stock et al., 2020).

2.4 Model spinup and simulation

Similar to Ross et al. (2023), we initialized the 1993–2019 hindcast simulation from rest starting at 1 January 1993, with ocean physics prescribed from GLORYS (described above), and we initialized the ocean biogeochemistry from a 10-year spinup simulation. We generated the spinup simulation by starting the model integration from rest on 1 January 1993 and by repeating ERA5 atmospheric conditions for 1993–1994 (May–December of 1993 and January–April 1994; following Stewart et al., 2020) for 10 1-year cycles. Atmospheric CO₂ was maintained as the 12-month, 1993 seasonal climatology, and the ocean boundaries were forced with a smoothed, daily climatology (i.e., averaged by “day of year” and smoothed with a triangular filter) of the hindcast's GLORYS12 1993–2019 open boundary conditions. River runoff was similarly prescribed as a smoothed daily climatology. The biogeochemical tracer fields at the end of this 10-year spinup simulation were then used to initialize biogeochemistry for the 27-year hindcast simulation.

The purpose of implementing a spinup was to omit drifts in the biogeochemistry associated with the adjustment of the model from its initialized state, which was generally based on coarse-resolution observation-based products, to the model's characteristic solution. We focused on fisheries-relevant variables in the top 500 m. We found that a spinup period of 10 years generally resolved initial model adjustments, which were strongest in the British Columbia region (Fig. S3 in the Supplement). While 10 years removed the strongest drifts, subtle trends remain in some regions, suggesting the potential value of longer spinup periods, particularly for representing the deeper ocean. These spinup sensitivities are left to future NEP10k development efforts.

2.5 Model evaluation

As described in Sect. 1, the model evaluation focuses on the simulation's capacity to represent fisheries and ecosystem-relevant features across and within the diverse ecosystems included within the NEP10k domain. The model evaluation therefore includes comparisons against both large-scale physical and biogeochemical patterns spanning the full domain (Sect. 2.5.1) and ecosystem-specific quantities (Sect. 2.5.2). These latter quantities were often drawn from Ecosystem Status Reports developed by NOAA fisheries to strategically inform marine resource management decisions (e.g., Ferriss, 2023; Siddon, 2023; Leising et al., 2024). Comparisons against spatial and seasonal patterns were complemented with interannual time series comparisons where possible; the latter serves as a building block toward making predictive applications. We note that several comparisons are made against gridded data products that were also used to force and initialize the NEP10k hindcast (i.e., GLORYS, TPXO, WOA23). While these comparisons are not fully independent, they are nonetheless meaningful tests of the capacity of the regional model to translate horizontal boundary and surface forcing into an interior solution that remains consistent with observations. The regional model must explain multiple observed interior properties by dynamically extending from the specified boundaries with a single set of self-consistent explicitly specified dynamics without the benefit of assimilating, or being informed by, observation from within the domain. Maintaining agreement with observation-based products in the domain interior thus supports the fidelity of these dynamics. We lastly assess the computational performance and viability of the model using analyses described in Sect. 2.5.3.

2.5.1 Full domain comparisons

We broadly evaluated NEP10k performance against gridded surface and 3D observation-based or observation-assimilated physical and biogeochemical products to assess the simulation's coastwide capacity to represent cross-ecosystem patterns. Table A2 summarizes these products and the time frames analyzed. For spatial comparisons and calculations, we first plot both the NEP10k results and the comparison product on their native grids using the Python geographic plotting package Cartopy (Met Office, 2022). We then regridded the finer-resolution product output (typically NEP10k, but not in the case of comparisons against GLORYS12 and chlorophyll comparisons) to the coarser-resolution comparison grid using the Earth System Modeling Framework (Hill et al., 2004) Python Regridding Interface (ESMPy) or xESMF conservative regridding (Zhuang et al., 2023). Unless otherwise stated, assessments include the area-weighted spatial mean bias (Bias, NEP10k – comparison data product), area-weighted root mean squared error (RMSE), median absolute error (MedAE), and Pearson correlation coefficient (R, based on the spatial pattern). We omit analysis of model performance in the Chukchi Sea (i.e., north of the Bering Strait at 66° N) – this region is included in the model integration due to the rectilinear nature of the grid and our objective to include the entire Bering Sea for which the Chukchi provides a boundary condition. However, it is not a primary region of interest for this model application and will be assessed in a nascent pan-Arctic MOM6 configuration.

For ocean temperature validations, we compared conditions against version 2.1 of the Daily Optimum Interpolation Sea Surface Temperature product (OISSTv2.1; Huang et al., 2021) and against GLORYS12 for both surface and subsurface conditions. OISSTv2.1 is generated from multiple temperature data sources and interpolated to a $\mathrm{1}/\mathrm{4}$° global grid, whereas GLORYS12 is a global eddying ($\mathrm{1}/\mathrm{12}$°) data-assimilative ocean reanalysis that demonstrates strong coherence with in situ surface and subsurface temperature records along the U.S. West Coast (Amaya et al., 2023a). Both reference products have a continuous monthly output covering 1993–2019.

NEP10k surface and subsurface salinity is compared against GLORYS12 reanalysis as well as the observation-based NOAA National Centers for Environmental Information (NCEI) $\mathrm{1}/\mathrm{10}$° Northern North Pacific (nnp; Version 2, Seidov et al., 2023) and Northeast Pacific (nep; Seidov et al., 2017) regional climatologies for salinity. Annual and seasonal means were downloaded for both nep and nnp regions for the decades 1995–2004 and 2005–2014 (the second decade for the older nep climatology only extends 2005–2012). To ensure temporal coherence, we regrid NEP10k separately for each region, using only the years represented by each regional climatology (i.e., 1995–2012 for the nep and 1995–2014 for the nnp). The two decadal, annual, and seasonal means for the regional climatologies are time-weight averaged, and then the regional climatologies and regridded NEP10k output are combined into a common grid. Where the nnp and nep regions overlap in the GOA (i.e., above 50° N), we use the values from the more recent nnp climatology.

We validated NEP10k mixed layer depth (MLD) against the 1° de Boyer Montégut (2024) monthly MLD climatology, which incorporates measurements from an assemblage of MBT, XBT, and CTD casts and profiling floats and defines the MLD as the seawater depth where the potential density is 0.03 (kg m³) greater than the density at a reference depth of 5 m. From NEP10k, we used the MOM6 diagnostic variable MLD_003, which calculates the mixed layer depth based on a user-defined reference depth (in our case, 5 m for consistency with de Boyer Montégut). The mixed layer depth is identified as the depth where the potential density increases by 0.03 kg m³ relative to the surface reference depth. We also compared NEP10k MLD against GLORYS12. The approximately equivalent MLD for GLORYS12 was determined by first calculating the potential density from monthly GLORYS12 potential temperature and salinity using the Python implementation of the Gibbs SeaWater (GSW) Oceanographic Toolbox of TEOS-10 (McDougall and Barker, 2011). We then calculated GLORYS12 MLD using the same criteria as de Boyer Montegut (2024) and the NEP10k MLD_003 diagnostic (i.e., depth at which potential density is 0.03 kg m³ greater than the density at 5 m depth at a given location).

NEP10k sea surface height (SSH) was compared against GLORYS12 sea surface height above geoid (zos) and absolute dynamic height (adt) above the Earth's geopotential surface (i.e., geoid) from 0.083° resolution satellite altimetry (CMEMS, 2023). Given the different reference frames for each observation, reanalysis, and model product, we mean-centered each dataset by subtracting its respective area-weighted time mean within the NEP10k region in order to facilitate direct comparison of seasonal and annual mean sea surface height distribution and gradients.

Tidal phase and amplitude for the M2 and K1 constituents were calculated using the hourly NEP10k sea surface height with the Unified Tidal Analysis and Prediction Python software package (Codiga, 2011). These tidal phases and amplitudes were compared against TPXO9 to demonstrate the ability of the model to incorporate and propagate tidal boundary forcings. We further included additional comparisons of tidal harmonics against several NOAA tide gauges (https://tidesandcurrents.noaa.gov/, last access: 8 May 2025) in the tidally complex eastern Bering Sea and western Gulf of Alaska.

NEP10k annual mean surface and subsurface nitrate and phosphate concentrations were compared against the 1° 2023 World Ocean Atlas (WOA23; Garcia et al., 2023a) for the time period 1993–2019. Primary phytoplankton nutrient limitation was calculated for annual and seasonal mean time frames following the methods detailed in Stock et al. (2020). These nutrient limitation distributions specifically illustrate where macronutrients nitrate and phosphate or micronutrient iron are the primary nutrient limitation of phytoplankton growth.

Surface chlorophyll is compared against the European Space Agency's satellite product produced as part of their Ocean Color Climate Change Initiative (OC-CCI; Sathyendranath et al., 2019, 2023). Monthly OC-CCI chlorophyll a fields from 1998 to 2019 were remapped from 4 km resolution to the coarser NEP10k grid. NEP10k grid cells where the OC-CCI satellite product is missing data were also masked in the corresponding month to ensure the annual and seasonal means are spatiotemporally consistent. Chlorophyll values were then log10 transformed before comparison.

We compared seasonal means of 200 m integrated mesozooplankton carbon biomass concentrations against the Coastal and Oceanic Plankton Ecology, Production and Observation Database (COPEPOD; Moriarty and O'Brien, 2013). As described in Ross et al. (2023), we scaled the COPEPOD dataset by a factor of 2 because the zooplankton represented in COBALT's mesozooplankton diagnostic (medium + large, ranging from 200 to 20 000 µm equivalent spherical diameter) likely represents a larger fraction of zooplankton biomass than in the COPEPOD observations, which are derived from collections that used a net mesh of 333 µm (Moriarty and O'Brien, 2013) and would exclude some of the size classes in the COBALT diagnostic (Skjoldal et al., 2013). This conversion is consistent with those typically found when comparing 200 µm and 333 µm mesh nets (Moriarty and O'Brien, 2013; Shropshire et al., 2020).

Similar to inorganic nutrients, surface and subsurface dissolved oxygen concentrations were compared against 1° WOA23 (García et al., 2023b) for 1993 through 2019, with NEP10k oxygen values being remapped to the WOA23 grid. We also computed the hypoxic boundary layer depth, here defined as the depth at which oxygen concentrations drop below 61.7 µmol O₂ per kilogram of seawater, as in Dussin et al. (2019).

We compared the annual and seasonal mean, surface and subsurface carbonate chemistry diagnostics, total alkalinity, dissolved inorganic carbon, and aragonite saturation state against corresponding values in the 1° Coastal Ocean Data Analysis Product in North America (CODAP-NA; Jiang et al., 2021) dataset (Jiang et al., 2022) for the period of 2004–2018.

2.5.2 Regional comparisons

The full domain comparisons were complemented with key fisheries-critical regional time series comparisons. While regions often have unique fisheries and ecosystem-critical patterns, temperature and chlorophyll variability are broadly important across ecosystems. We thus complemented the broad spatial comparisons with region-specific time series of shelf (defined as grid cells where the bottom depth is less than 500 m) conditions, where the subregions are those shown in Fig. 1 and regional shelf extents are depicted in Fig. S2 in the Supplement. Both monthly climatologies and anomaly (with the 12-monthly climatological cycle removed) time series for surface and bottom temperatures were compared against GLORYS12, while time series of chlorophyll were compared against OC-CCI. For these (and later) time series analyses, we report the Pearson correlation coefficient within the respective figure as well as the Kling–Gupta efficiency (KGE; Gupta et al., 2009) and its components in the Supplement (Table S1) for a more comprehensive assessment of the interactions of time series correlation, bias, and variance. It should be noted that the KGE is calculated using the full time series rather than the climatology or the anomaly time series and thus the Pearson correlation coefficients may differ between the figures and the supplemental table.

For additional environmental context, anomaly time series are depicted against warm and cold episodes of the Ocean Niño Index published by the NOAA Climate Prediction Center (https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php, NOAA Climate Prediction Center, 2023a), where the warm and cold episodes are defined as periods when the 3-month running mean of sea surface temperature (SST) anomaly in the Niño3.4 region is above or below 0.5 °C, respectively. The purpose of this comparison is to ascertain whether the model is able to accurately recreate the strength of the relationship between local variability and this foremost mode of global climate variability. Variations in simulation skills for different depth ranges within each subregion were also analyzed to assess changes in model fidelity in more inshore and offshore regions.

Additional region-specific assessments are described for the Bering Sea, Gulf of Alaska, and California Current below. Given the length constraints of a single documentation paper, we limited treatment to two to three of the most prominent ecosystem indicators currently used for each system beyond the foundational temperature and chlorophyll comparisons described above.

Our additional evaluation in the Bering Sea focused on the representation of the Bering Sea cold pool and sea ice extent. As discussed in Sect. 1, fluctuations in the bottom area covered by the Bering Sea cold pool, generally defined as waters with <2 °C in the summer (Wyllie-Echeverria and Wooster, 1998; Mueter and Litzow, 2008), have been associated with a range of ecosystem impacts (e.g., Clement Kinney et al., 2022). Cold pool dynamics are intertwined with sea ice fluctuations, with sea ice also having important implications for the timing of seasonal ecosystem transitions (Wyllie-Echeverria and Wooster, 1998; Mueter and Litzow, 2008; Brown and Arrigo, 2013; Hunt et al., 2022).

For the Bering Sea cold pool, we spatially and temporally interpolated the daily NEP10k bottom temperature using the Python package xESMF (Zhuang et al., 2023) to correspond with Alaska Fisheries Science Center (AFSC) Bottom Trawl Survey gear temperature samples collected from 1993 to 2019. These data are available in the Alaska Fisheries Science Center cold pool GitHub repository (https://github.com/afsc-gap-products/coldpool, NOAA-AFSC, 2024a). We compared the trawl survey station bottom temperatures from the NEP10k simulation against the AFSC dataset following the methods in Kearney (2021) and analyzed the interpolated model output using the cold pool toolset to reproduce the cold pool area (CPA) indices reported by Rohan et al. (2022).

We compared seasonal Bering Sea sea ice against satellite observations from the National Snow and Ice Data Center (NSIDC; dataset NSIDC0051; Cavalieri et al., 1996). We also compared both the spatial mean extent in the entire Bering Sea and temporal coherence in the southeastern Bering Sea.

Hauri et al. (2024) highlight how the interaction of different localized modes of multiannual to decadal climate variability can predispose the Gulf of Alaska to extreme physical and biogeochemical events. These climate variations are most visibly reflected in observed Gulf of Alaska SSH variability. The first principal component of the detrended and deseasonalized SSH over the Gulf of Alaska (62° N 50° N, 160° W 135° W) was referred to as the Northern Gulf of Alaska Oscillation (NGAO, Hauri et al., 2021). A positive phase is associated with weak cyclonic winds over the subpolar gyre, resulting in a higher SSH and decreased Ekman-driven upwelling (i.e., Ekman suction). This state is associated with warmer temperatures but reduced prevalence of deep high-acidity water. That is, risks of thermal stress are enhanced, while risks of acidification stress are reduced, with the opposite effects for negative NGAO. The second principal component of the detrended and deseasonalized SSH variability is referred to as the Gulf of Alaska downwelling index (GOADI; Hauri et al., 2024). The GOADI serves as a proxy of downwelling strength for Gulf of Alaska coastal waters: a positive index is associated with elevated coastal SSH, enhanced coastal downwelling, and a reduced risk of the intrusion of cold, acidic, and low-oxygen water onto the bottom of the Gulf of Alaska shelf. This intrusion risk is heightened under negative GOADI.

We assessed NEP10k's ability to generate realistic NGAO and GOADI patterns by comparing against satellite altimetry from the Copernicus Marine Environment Monitoring Service (CMEMS, 2023). Empirical orthogonal function analysis was performed on SSH across the GOA domain in a manner consistent with Hauri et al. (2021) and Hauri et al. (2024). We then generated composites of ecosystem conditions during the positive vs. negative phases of the GOADI to assess whether NEP10k can successfully recreate the shelf-scale surface and benthic condition anomalies that significantly impact living marine resource habitat and well-being (Hauri et al., 2024).

Fisheries and ecosystems in the California Current are shaped by the timing, strength, and source waters fueling the strong seasonal upwelling. The system-specific indicators chosen for this region thus focus on these patterns. First, we compared the vertical mass transport (calculated as the depth-integrated divergence of orthogonal horizontal mass transports) at a depth of 30 m to the Coastal Upwelling Transport Index (CUTI) developed by Jacox et al. (2018). As in Jacox et al. (2018), transports were integrated to 75 km offshore over 1° latitude bins. We assessed long-term trends in dissolved oxygen concentrations against those calculated at stations in the California Cooperative Oceanic Fisheries Investigations (CalCOFI) observation array similar to the methods of Bograd et al. (2008). We interpolated monthly 3D NEP10k dissolved oxygen to the locations and depths of the CalCOFI bottle sample data (https://calcofi.org/data/oceanographic-data/bottle-database/, last access: 7 February 2025) from 1993 to 2019. We then calculated linear trends for both NEP10k and CalCOFI at specific station locations. We also included additional comparisons of NEP10k representation of CalCOFI temperature, salinity, and biogeochemistry measurements.

2.5.3 Computational expense and scaling

As mentioned in Sect. 2.2, simulations were conducted on NOAA's GAEA high-performance computing system. This system consists of HPE-Cray EX 3000 nodes (2 × AMD EPYC 9654, 2.4 GHz base, 96 cores per socket), connected via HPE Slingshot 11 – a high-speed interconnect designed for exascale systems. The system also features over 150 PB of shared storage using IBM Spectrum Scale parallel file systems. The model runs in a distributed-memory configuration using MPI across hundreds to thousands of cores. Additional system details can be found in the NOAA RDHPCS documentation (https://docs.rdhpcs.noaa.gov/systems/gaea_user_guide.html#system-overview, last access: 22 October 2024).

As described in Sect. 1, the viability of the NEP10k configuration for ecosystem applications depends on its ability to not only simulate fisheries-critical features but also to run with sufficient computational economy to permit generation of the thousands of years of retrospective forecasts and projections required to provide credible uncertainty estimates (e.g., Koul et al., 2024; Ross et al., 2024). However, we also recognize that others interested in running the NEP10k configuration may have different computing resource availability. Therefore, we report the computational performance under different NEP10k configuration options (i.e., scaling, land masking, and time-step splitting) in order to provide insight into how one might optimize production on a given computing system.

To quantify computational performance, we focused on the scaling of the wall clock time for 1 year of simulation against the number of processing elements (PEs). Variations in both the number and layout of PEs were considered. For our baseline production simulations herein, we divided the NEP10k domain (342 columns × 816 rows of tracer grid cells) across 32×80 PEs. This division yields an $\sim \mathrm{10}\times \mathrm{10}$ grid (i.e., square) decomposition of model grid cells on each PE. Land processor masking in MOM6 further economizes computational resources by omitting domain subregions without ocean (i.e., those that contain only land) grid cells from PE assignment, thus presenting a domain-specific optimization consideration when selecting a specific PE configuration. We were able to mask 524 PEs with the 32×80 PE breakdown, so our total PE count for this configuration was 2036 (20 % fewer than the otherwise 2560 PEs required for this breakdown).

The scalability of the simulation with increasing and decreasing processor counts was explored using alternative layouts with fewer PEs (40×40), a similar PE total but with a more rectangular model grid cell decomposition (a 50×50 PE breakdown yielding an $\sim \mathrm{7}\times \mathrm{16}$ model grid cell subset per PE), and larger numbers of PEs (50×75 and 50×100). These experiments allow us to judge the relative efficiency of our base configuration and the point of diminishing returns as the PE count is increased and growing requirements for inter-PE communication begin to overwhelm the advantage of more PEs. Finally, we include additional 50×75 PE and 50×100 PE simulations with the thermodynamic time step equal to the baroclinic time step (400 s) rather than 3 times the baroclinic time step (i.e., 1200 s), as was used in the base configuration. These last two experiments allow us to quantify and demonstrate the computational value of the flexible time stepping that MOM6 enables.

3 Results

3.1 Domain-wide evaluation

3.1.1 Large-scale physical ocean properties

Annual mean SST and subsurface temperatures broadly agree with the distribution and curvature of reference isotherms along the west coasts of the United States, Canada, and Mexico (Fig. 2), with temperatures largely falling within 0.5 °C of OISST (Fig. 2c, RMSE=0.28 °C) and GLORYS12 SST values (Fig. 2f, RMSE=0.29 °C). A surface temperature cold bias of just over 0.5 °C is apparent over the eastern Bering Sea, while a warm bias of similar magnitude is apparent in the nearshore regions of the southern and central California Current System. At a depth of 200 m, larger warm biases relative to GLORYS12 are apparent in the Gulf of Alaska, where the northern edge of the eastward flowing North Pacific Current interacts with the adjacent westward flowing Alaska Stream (Fig. 2l, Stabeno et al., 2004), and a warm bias of similar magnitude appears in the southwest corner of the domain. These biases are seasonally persistent during both Boreal winter (January–March, Fig. S1 in the Supplement) and summer (July–September, Fig. S2), as are the cold (Fig. S1c and f) and warm (Fig. S2c and f) coastal surface biases, respectively. In all seasons and across depths above 200 m, however, the overall absolute model bias is below 0.38 °C, the RMSE stays below 0.57 °C, and the correlations with OISSTv2.1 and GLORYS12 stay above 0.98 (Figs. 2, S1, and S2).

Figure 2 Temperature comparisons. Annual mean surface and subsurface (100 m, 200 m) temperature compared against NOAA OISSTv2.1 and the GLORYS12 reanalysis. Values in the left two columns represent the average of the annual means covering 1993 through 2019. The right column depicts the difference between NEP10k and the respective validation product, along with the area-weighted mean bias, root mean squared error (RMSE), medium absolute error (MedAE), and Pearson correlation coefficient (R). The NEP10k model domain below 66° N is outlined in black. Panels a and d show the same model output.

Similar to temperature, NEP10k broadly reproduces annual mean salinity fields found in regional climatologies and GLORYS12, with the majority of the domain falling within 0.25 practical salinity units (PSU) of the reference datasets (Fig. 3). Notable fresh surface biases exceeding 0.5 PSU occur along the coast in the Gulf of Alaska, eastern Bering Sea, and northern CCS, coincident with regions of substantial freshwater inputs from rivers and glacial melt (Fig. 3c and l). Positive salinity biases relative to GLORYS12 occur in the western Bering Sea at the surface and 100 m and over all depths in the southwest region of the domain (Fig. 3, right panels). In the latter case, the salty bias coincides with warm biases (Fig. 2). Seasonally, similar generally modest biases can be seen in the Boreal winter (Fig. S3) and summer (Fig. S4 in the Supplement) equivalents.

Figure 3 Salinity comparisons. Annual mean surface and subsurface (100 m, 200 m) salinity compared against NCEI regional ocean climatologies and the GLORYS12 reanalysis. The regional climatologies are a composite of the northeast Pacific (nep) and northern North Pacific (nnp) climatologies. The nep climatology extends from 1995 to 2012, while the updated nnp climatology (Version 2) covers 1995–2014. Where the two regional climatologies overlap in the GOA (i.e., above 50° N), we use the more recent nnp climatology. For comparison against the model, we use the same years of NEP10k, with panels (a), (d), and (g) showing the model values for average annual mean salinities for 1995–2014 above 50° N (as opposed to average annual mean salinities for 1995–2012 below 50° N). The comparison against GLORYS12 (bottom three rows) covers 1993–2019. Area-weighted bias, root mean squared error (RMSE), median absolute error (MedAE), and Pearson correlation coefficient (R) are reported in the right column of figures, depicting the difference between NEP10k and the respective validation product.

Mixed layer depth in NEP10k, defined as the depth at which density is 0.03 kg m⁻³ greater than at 5 m depth, exhibits a modest shallow/negative bias relative to the estimates of de Boyer Montégut (2024), with deeper (positive) biases occurring in the interior ocean near the Bering shelf break (Fig. 4, top row). These biases are amplified and reduced during Boreal winter (JFM, Fig. S5 in the Supplement, top row) and summer (JAS, Fig. S6, top row), respectively, when mixing drivers (i.e., surface heating/cooling, wind, and storm intensity) are correspondingly modified. Conversely, NEP10k exhibits a positive mean bias when compared against GLORYS12 MLD, which is particularly pronounced in the Bering Sea (Fig. 4, bottom row) and exhibits a reverse seasonal response (i.e., reduced positive bias in the winter and increased in the summer; Figs. S5 and S6 in the Supplement, bottom row). With the exception of the deep/positive winter biases in the Bering Sea, the model represents MLD spatial variability fairly well, with significant (p<0.001) correlations exceeding 0.85 across all seasons and comparisons (Figs. 4, S5, and S6).

Figure 4 Mixed layer depth comparisons. Climatological mean of mixed layer depth compared against de Boyer Montégut (a–c) and GLORYS (d–f). Black reference contours in (a), (b), (d), and (f) are depicted at 5 m intervals and at 8 m intervals in (c) and (f); contours depicting negative values in (c) and (f) are drawn with dashed lines. Area-weighted bias, root mean squared error (RMSE), median absolute error (MedAE), and Pearson correlation coefficient (R) are reported in the right column of figures, depicting NEP10k-respective reference products. All values represent the annual mean for years 1993 through 2019, and the extent of the NEP10k domain is outlined in black in all figures. Panels a and d show the same model output.

SSH gradients in the NEP10k hindcast are broadly consistent with GLORYS12 and CMEMS satellite altimetry (Fig. 5), exhibiting the lowest values along the Aleutian island chain, in the GOA, and in the western Bering Sea and the highest values near 25° N along the western edge of the domain. Similar to satellite measurement and GLORYS12, NEP10k also exhibits relatively low SSH along the U.S. West Coast (compared with offshore SSH values at the same latitude), a signature of coastal upwelling. However, the SSH gradients in NEP10k are smaller along the Aleutian island chain than exhibited in the reference datasets. There is a notable correspondence of this SSH gradient bias with the Gulf of Alaska subsurface temperature biases noted in Fig. 2, suggesting a potential relationship between these two features.

Figure 5 Sea surface height comparisons. NEP10k average-centered, climatological mean sea surface height comparison for NEP10k (a, d; identical panels), GLORYS12 (b), CMEMS satellite altimetry (e), and their respective differences (c, f). All values represent the annual mean (1993–2019). Area-weighted mean bias (Bias), root mean squared error (RMSE), median absolute error (MedAE), and Pearson correlation coefficient (R) are reported in the right column of figures, depicting the difference between NEP10k and the comparison product; all correlations are significant (p<0.001). The reference height contours in all panels are drawn at 0.1 and 0.05 m intervals for the mean and difference plots, respectively, with negative values shown as dashed lines. All panels show the extent of the NEP10k domain as a black outline.

Compared against the TPXO dataset, which was used as the tidal boundary forcing conditions, NEP10k reproduces tidal amplitude and phases in the domain interior with high fidelity (Fig. 6). The greatest tidal amplitude discrepancies occur in the nearshore regions of the eastern Bering Sea (Fig. 6c and f) and partially enclosed features (e.g., northern Gulf of California and Cook inlet; Fig. 6c). Amplitude biases for the most prominent semidiurnal (M2) and diurnal (K1) constituents in these nearshore and partially enclosed regions can exceed 20 cm and 10 cm, respectively. These regions, however, also have the largest overall amplitudes, with values exceeding 1 m and 50 cm, respectively. Such nearshore tidal biases are not surprising given the relatively coarse 10 km resolution enlisted herein, and we note that skillful tidal simulations extend all the way to the coast in most regions. To investigate some of these biases further, we include additional zoomed-in maps of the eastern Bering Sea and western Gulf of Alaska in the Supplement (Fig. S12), along with comparison against several tide gauges in that region. Both TPXO and NEP10k perform well at most tide gauges. Generally, TPXO better approximates tidal harmonic constituents than NEP10k, with higher Pearson correlation coefficients and/or lower RMSE (with the exception of the M2 phase). However, in cases such as the gauge in Anchorage, AK, the bias in M2 amplitude for TPXO is comparable to the bias exhibited by NEP10k. Because these biases have opposite signs, the discrepancy between the two gridded products (i.e., NEP10k-TPXO, shown in the maps in Figs. 6 and S12) exaggerates the model bias by almost a factor of 2 relative to the bias for the gauge. Thus, some of the more severe nearshore differences in Fig. 6 may be a reflection of how NEP10k and TPXO approximate complex coastline geometry (bottom of Fig. S13 in the Supplement) rather than an exact indication of NEP10k performance.

Figure 6 M2 and K1 tidal amplitudes and period. Comparison of tidal constituents M2 (top row) and K1 (bottom row) in NEP10k against those in the TPXO9 forcing dataset. Filled contours depict the tidal amplitude, while overlain colored contours depict the tidal phase for the given constituent. Filled contours in the difference plot (c, f) show the difference in amplitude only; bias, root mean squared error (RMSE), median absolute error (MedAE), and Pearson correlation coefficient (R) are also reported in these panels. The extent of the NEP10k domain is outlined in gray in all figures.

3.1.2 Large-scale biogeochemical and ecosystem properties

Macronutrient concentrations (nitrate and phosphate) exhibit large-scale agreement with annual World Ocean Atlas nutrients, but significant regional biases are also apparent (Figs. 7 and 8). The largest high bias occurs along the Aleutian island chain and Bering Sea shelf break. In the simulation, the region of elevated surface nutrients observed in the central Bering Sea extends further south and east in the model. These biases correspond with the most prominent region of overmixing (Fig. 4). Positive surface nitrate and phosphate biases in affected regions exceed 5 µmol kg⁻¹ NO₃ and 0.25 µmol kg⁻¹ PO₄, respectively, and extend with lesser severity onto the Bering shelf. The positive surface bias is underlain by negative nitrate and phosphate biases at 200 m, reinforcing the likelihood that the surface high macronutrient bias is linked to excessive mixing rather than excessive nutrients in underlying source waters. Uncertainty in nitrogen removal processes in shallow Bering shelf sediments (e.g., denitrification and burial) may also play a role in the perpetuation of biases onto the shelf. Macronutrient concentrations in Gulf of Alaska surface waters, in contrast, are biased low by 1.5–3 µmol kg⁻¹ NO₃ and 0–0.375 µmol kg⁻¹ PO₄, respectively (Figs. 7c and 8c), despite exhibiting a combination of positive and negative biases at depth. These biases are consistent with shallow mixed layer biases in the Gulf of Alaska (Fig. 4). Finally, the California Current exhibits a modest positive surface macronutrient bias. Despite these discrepancies, the simulation generally exhibits high correlations with observed macronutrients (R>0.96) and RMSEs that are only ∼5 % of the dynamic range of the macronutrient concentrations across the west coast ecosystems. This skill extends to seasonal patterns with correlation values exceeding 0.8 and RMSE<10 % of the dynamic range in all cases (Figs. S9–S12 in the Supplement). Notably, winter and summer nitrate conditions exhibit more pronounced bias patterns relative to the mean state, with particularly high levels in the Bering surface waters and low levels in portions of the Gulf of Alaska (Figs. S9c and S10c). Conversely, surface phosphate levels over the Bering shelf are biased low in the winter and high in the summer (Figs. S11c and S12c). Summer surface nitrate levels along the California Current Ecosystem (CCE) (Fig. S10c) are potentially suggestive of overrepresentation of summer upwelling.

Figure 7 Nitrate comparisons. Annual mean surface and subsurface (100 m, 200 m) nitrate compared against WOA23. Comparison time frames cover 1993–2019. Reference contours are depicted in black at 5 and 1.5 µmol nitrate per kg sea water in the mean state (left and center columns a, b, d, e, g, h) and difference (right column c, f, i) plots; contours representing negative values in the difference plot are drawn as dashed lines. Bias, root mean squared error (RMSE), median absolute error (MedAE), and Pearson correlation coefficient (R) are reported in the right column (c, f, i) of figures, depicting the difference between NEP10k and WOA23. The extent of the NEP10k domain is outlined in black in all figures.

Figure 8 Phosphate comparisons. Annual mean surface and subsurface (100 m, 200 m) phosphate compared against WOA23. Comparison time frames cover 1993–2019. Reference contours are depicted in black at 0.25 µmol phosphate per kg sea water in the mean state (left and center columns a, b, d, e, g, h) and difference (right column c, f, i) plots; contours representing negative values in the difference plot are drawn as dashed lines. Bias, root mean squared error (RMSE), median absolute error (MedAE), and Pearson correlation coefficient (R) are reported in the right column (c, f, i) of figures, depicting the difference between NEP10k and WOA23. The extent of the NEP10k domain is outlined in black in all figures.

While macronutrients play an important role in the biogeochemistry and ecosystem dynamics of the NEP, iron has been observed to be a limiting or co-limiting nutrient (Browning et al., 2017; Browning and Moore, 2023). The simulated distribution of surface iron exhibits a gradient from inshore highs exceeding 1 nmol kg⁻¹ to offshore lows <0.25 nmol kg⁻¹ (Fig. 9, left panel). This distribution of dissolved iron results in large-scale patterns of phytoplankton iron limitation in the NEP10k simulation (Fig. 9, right panel) that are consistent with those observed (e.g., Moore et al., 2013; Hutchins et al., 1998).

Figure 9 Surface dissolved iron and phytoplankton nutrient limitation. NEP10k simulated annual mean surface dissolved iron concentrations (a) and climatological mean distribution of the nutrient most limiting to phytoplankton growth (b). In COBALT, the degree of limitation by N, P, and Fe is expressed as a factor between 0 and 1 (Stock et al., 2020). Nutrient limitation is then calculated according to Liebig's law of the minimum. This most limiting nutrient is indicated in the figure below. We further differentiate areas where the N, P, or Fe limitation term is less than 0.25 more limiting than another nutrient, which effectively indicates areas that are near co-limitation. Time frame covers 1993–2019. Note: Sparse P limitation occurs nearshore.

Simulated surface chlorophyll is spatially well correlated with satellite-based chlorophyll estimated from the OC-CCI (Fig. 10), and simulated values are generally within a factor of 2 of those observed, which span 2 orders of magnitude (i.e., the RMSE of the log₁₀-transformed data is less than 0.3 in all seasons). The simulation, however, is generally biased high in the Gulf of Alaska and Bering Sea in the boreal spring and summer, with biases exceeding a factor of 2 along the Bering Sea shelf break and along the subpolar/subtropical boundary in the Gulf of Alaska. The model underestimates the OC-CCI-based chlorophyll concentration during the fall and winter on the eastern Bering Sea shelf: while NEP10k-COBALTv3 suggests lower chlorophyll concentrations during these cold and dark periods, OC-CCI estimates remain high in nearshore waters. Indeed, satellite-based estimates suggest higher chlorophyll along the Bering coast in fall and winter than in spring and summer. It is notable, however, that satellite-based chlorophyll estimates are sporadic at high latitudes during these seasons, and OC-CCI uses a chlorophyll estimation algorithm developed primarily for “case 1”/oceanic water. Vigorously mixed, turbid waters along the Bering shelf in winter undoubtedly depart considerably from the algorithm's high degree of water transparency assumptions. In the CCE, the model is able to match the juxtaposition of coastal chlorophyll highs and subtropical offshore lows estimated by OC-CCI during the spring and summer upwelling period. Elevated chlorophyll levels do extend further offshore in the simulation than satellite estimates suggest. Values are also elevated near the domain boundary during this period, likely due to some spurious boundary mixing. Fall and winter conditions in the California Current exhibit a moderate positive bias in offshore waters that generally falls below a factor of 2.

Figure 10 Surface chlorophyll comparisons. Seasonal means of surface chlorophyll compared with OC-CCI satellite observations. The 3-month seasonal periods include January through March (JFM, a–c), April through June (AMJ, d–f), July through September (JAS, g–i), and October through December (OND, j–l). Comparison time frames cover 1998–2019; all chlorophyll values were log₁₀-transformed prior to temporal averaging. Bias, root mean squared error (RMSE), median absolute error (MedAE), and Pearson correlation coefficient (R) are reported in the right column (c, f, i, l) of figures, depicting the difference between NEP10k and OC-CCI. Black contours in the right column (c, f, i, l) indicate where the difference $=\pm {\log }_{\mathrm{10}}\left(\mathrm{2}\right)$. The extent of the NEP10k domain is outlined in gray in all figures.

Moving up the food web, simulated seasonal mesozooplankton biomass concentrations (Fig. 11) exhibit similar large-scale spatial and seasonal patterns as the COPEPOD database (Moriarty and O'Brien, 2013). The patchiness of the observations reduces correlations relative to the smoother physical, nutrient, and satellite-based chlorophyll estimates compared thus far (R≥0.30 for all seasons). However, peak summer concentrations (∼50 mg C m⁻³) consistent with observed values are evident in the Bering Sea and inshore regions of the Gulf of Alaska in both the model and observations. These highs contrast sharply with observed and modeled values (∼1–2 mg C m⁻³) within the North Pacific subtropical gyre. Intermediate values of ∼10–20 are evident in the California Current upwelling. Both the observed and modeled values are highest during the peak summer upwelling period, though the highest modeled values are somewhat lower, particularly in nearshore regions. This pattern will be addressed further in the Discussion section. The offshore waters of the Gulf of Alaska and western Bering Sea exhibit summer mesozooplankton biomass peaks of similar magnitude as the California Current, with simulated values again lower yet comparable to those observed.

Figure 11 Seasonal zooplankton biomass. Seasonal mean mesozooplankton biomass concentrations for NEP10k on the model grid (top row), the COPEPOD dataset (middle row), and NEP10k values remapped to the COPEPOD grid where there are corresponding data from the COPEPOD dataset (bottom row). The bottom row also reports statistics using the log₁₀ normalized data, specifically the area-weighted mean bias (Bias, NEP10k – COPEPOD), the area-weighted root mean squared error (RMSE), the median absolute error (MedAE), and the Pearson correlation coefficient (R); all correlation values are significant (p<0.001). Maps are plotted with a gray background to increase contrast with the patchy observation data.

Simulated oxygen concentrations in the top 200 m in the NEP10k are generally spatially consistent with WOA (Fig. 12). Some biases, however, are apparent below the surface. Most notably, the model has a low oxygen bias south of the Aleutian Islands at 100 m (Fig. 12f). This bias coincides with a warm water bias (Fig. 2) and is overlain by a fresh/high stratification bias (Figs. 3 and 4). As noted above, this is the region where the westward flowing Alaska Stream and eastward flowing North Pacific Current interact, suggesting that the biases may be linked to a suboptimal representation of these two currents. Moderately high oxygen biases (i.e., greater than 25 µmol kg⁻¹) are apparent in the western Bering Sea, eastern Gulf of Alaska, and off of Baja at 200 m (Fig. 12i), but none are large enough to compromise NEP10k's large-scale fidelity to the observed oxygen distribution in the top 200 m (i.e., R values≥0.9 across depths and seasons; Figs. 12, S14, and S15 in the Supplement).

Figure 12 Dissolved oxygen comparisons. Annual mean surface and subsurface (100 m, 200 m) dissolved oxygen compared against WOA23. Comparison time frames cover 1993–2019. Reference contours are depicted in black at 25 µmol oxygen per kg sea water in the mean state (left and center columns a, b, d, e, g, h) and difference (right column c, f, i) plots; contours representing negative values in the difference plot are drawn as dashed lines. Bias, root mean squared error (RMSE), median absolute error (MedAE), and Pearson correlation coefficient (R) are reported in the right column (c, f, i) of figures, depicting the difference between NEP10k and WOA23. The extent of the NEP10k domain is outlined in black in all figures.

Deeper in the water column, NEP10k robustly simulates the cross-ecosystem variation in the depth of the hypoxic boundary (i.e., the depth at which the oxygen concentration drops below 61.7 µmol oxygen per kg sea water; see Fig. 13). The hypoxic boundary is shallowest, approaching 100 m from the surface, along the southern domain boundary, which lies along the periphery of the broader eastern equatorial Pacific hypoxic zone. The hypoxic boundary then descends progressively to ∼400 m in both the model and observations as one moves northward along the California Coast into Canada, before shoaling again to ∼150 m in the northern Gulf of Alaska. While these overall patterns are consistent, the biases discussed in Fig. 12 are echoed in the hypoxic boundary layer depth. The boundary layer is deeper in the western Bering Sea, eastern Gulf of Alaska, and Southern CCS but biased shallow south of the Aleutian island Chain and, to a lesser degree, in the northern-to-central CCS.

Figure 13 Hypoxic boundary layer depth. Annual mean hypoxic boundary layer depth (i.e., depth at which the dissolved oxygen concentration drops below 61.7 µmol oxygen per kg sea water) compared against WOA23. Black reference contours indicate 150 m and 25 m intervals in the mean state (a, b) and difference (c) plots; contours representing negative values in (c) are drawn as dashed lines. Area-weighted mean bias (Bias), root mean squared error (RMSE), median absolute error (MedAE), and Pearson correlation coefficient (R) are reported in panel (c). The extent of the NEP10k domain is outlined in black in all figures.

Finally, simulated carbon chemistry patterns (total alkalinity, dissolved inorganic carbon (DIC), and aragonite saturation state; Figs. 14–16) broadly capture observation-based estimates reported in CODAP-NA. Low coastal surface alkalinity patterns consistent with low alkalinity river inputs are apparent in the Gulf of Alaska and, to a lesser degree, the eastern Bering Sea. Simulated alkalinity increases from these lows toward maximal values in the North Pacific gyre in a manner consistent with observations, though the simulated values are biased high (Fig. 14a–c). The largest positive surface alkalinity biases occur in the western Bering Sea and in the southwest corner of the domain. These surface alkalinity biases are aligned with positive salinity biases that penetrate to depth (Fig. 3). The largest subsurface bias, however, occurs at a depth of 100 m in the Gulf of Alaska near the large freshwater outflows in the Gulf of Alaska. This bias distribution suggests that the low alkalinity freshwater signal in this region may be overly restricted to the surface in the model, though there does not appear to be a strong positive subsurface salinity model bias in this region (Fig. 3).

Figure 14 Total alkalinity comparisons. Annual mean surface and subsurface (100 m, 200 m) total alkalinity compared against CODAP-NA. Comparison time frames cover 2004–2018. Reference contours are depicted in black at 25 µmol alkalinity per kg sea water in the mean state (left and center columns a, b, d, e, g, h) and difference (right column c, f, i) plots. Area-weighted mean bias (Bias), root mean squared error (RMSE), median absolute error (MedAE), and Pearson correlation coefficient (R) are reported in the right column (c, f, i) of the difference plots. All correlation values are significant at p<0.001. The extent of the NEP10k domain is outlined in black in all figures.

Figure 15 Dissolved inorganic carbon comparisons. Annual mean surface and subsurface (100 m, 200 m) concentration of dissolved inorganic carbon compared against CODAP-NA. Comparison time frames cover 2004–2018. Reference contours are depicted in black at 50 and 25 µmol carbon per kg sea water in the mean state (left and center columns a, b, d, e, g, h) and difference (right column c, f, i) plots; contours representing negative values in the difference plots are drawn as dashed lines. Area-weighted mean bias (Bias), root mean squared error (RMSE), median absolute error (MedAE), and Pearson correlation coefficient (R) are reported in the right column (c, f, i) of the difference plots. All correlation values are significant at p<0.001. The extent of the NEP10k domain is outlined in black in all figures.

Figure 16 Aragonite saturation state comparisons. Annual mean surface and subsurface (100 m, 200 m) aragonite saturation state compared against CODAP-NA. Comparison time frames cover 2004–2018. Reference contours are depicted in black at 0.5 and 0.25 saturation state units in the mean state (left and center columns a, b, d, e, g, h) and difference (right column c, f, i) plots; contours representing negative values in the difference plots are drawn as dashed lines. Area-weighted mean bias (Bias), root mean squared error (RMSE), median absolute error (MedAE), and Pearson correlation coefficient (R) are reported in the right column (c, f, i) of the difference plots. All correlation values are significant at p<0.001. The extent of the NEP10k domain is outlined in black in all figures.

Dissolved inorganic carbon has a high bias that is consistent with the high alkalinity bias (compare Figs. 14 and 15). Like alkalinity, the largest positive biases occurred along the Bering Sea shelf break and in the southwestern corner of the domain where areas are overmixed (Fig. 4) and exhibit salty biases (Fig. 3). The high surface DIC bias in the northern Gulf of Alaska, however, is more pronounced than the corresponding high surface alkalinity bias in this region (i.e., Fig. 13c vs. Fig. 14c). The northern Gulf of Alaska is strongly impacted by river and glacial outflows. While some of these freshwater sources (e.g., the Copper and Susitna rivers) have observational constraints on DIC and Alk, most do not. Improved constraints may be needed to improve the model fit in this region.

The more-pronounced high surface DIC bias in the northern Gulf of Alaska yields aragonite saturation states that are 0.25–0.5 units lower than the CODAP-NA product (Fig. 16). The overall gradient between low saturation states (higher acidification vulnerability) in the surface waters of the Bering Sea/Gulf of Alaska to high saturation states (lower acidification vulnerability) in equatorial and subtropical surface waters in the southern parts of the domain, however, is well captured (Fig. 15c, R=0.93). Saturation state biases are also small in subsurface waters where subsaturated waters are more prevalent (Fig. 16, middle and bottom panel) and where valuable shell, crab, and demersal fisheries reside.

3.2 Region-specific evaluation

The evaluation of NEP10k against observed large-scale physical and biogeochemical patterns in Sect. 3.1 was generally favorable. In all cases, the model was able to capture the primary physical, biogeochemical, and plankton contrasts across ecosystems within the broad NEP10k domain with often high but at least moderate fidelity. As described in Sect. 1, however, the NEP10k configuration is intended for marine resource applications both across and within NEP10k subregions and across management-relevant time horizons from seasons to multiple decades. The evaluation in Sect. 3.1 provides a foundation for such applications but is not sufficient. The evaluation in this section focuses on regional fisheries-critical metrics and their variation across management-relevant seasonal to multidecadal time horizons.

Perhaps the most ubiquitous indicators of ecosystem state across all regions are ocean temperature (surface and bottom) and surface chlorophyll. These indicators are highly relevant to diverse aspects of ecosystem function, and long time series of observation-informed estimates are available. Modeled shelf (where depth <500 m) surface and bottom temperature climatologies for the regions identified in Fig. 1 exhibit high correlation (Fig. 17, left column) with GLORYS12, but surface temperatures tend to be biased warm in more southerly regions. As initially illustrated in Figs. 2 and S2, mean and summer surface temperatures, respectively, in the central and southern California Current System are 1–2 °C warmer than those observed, but biases in other regions tend to be <1 °C.

Figure 17 Surface and bottom temperature comparisons for shelf (0–500 m) regions. Regional shelf (depth ≤500 m) surface and bottom temperature climatologies (left column) and anomaly time series (right column) for the subregions delineated in Fig. 1. Comparison of temperature climatologies (left panels) and monthly anomalies (right panels) for surface (orange) and bottom (purple) temperatures for NEP10k (bold) and GLORYS12 (pale). Axes for surface and bottom temperature anomalies are separate and offset for improved readability. Pearson correlation coefficients are reported for surface (R_STC, R_STA) and bottom (R_BTC, R_BTA) climatology and anomaly comparisons, respectively. Background shading in the monthly anomaly time series plots indicates the oceanic Niño index produced by the NOAA Climate Prediction Center for context.

The NEP10k and GLORYS12 monthly surface and bottom temperature anomaly time series (Fig. 17, right column) have correlations >0.7 in nearly all regions, with values exceeding 0.9 in many. In the California Current, fluctuations in both NEP10k and GLORYS12 show a strong correspondence with the Niño 3.4 index (shaded regions), with warm conditions prevalent during warm ONI states and cold conditions prevalent during cold ONI. The lowest NEP10k-GLORYS12 correlations (R=0.82 for the surface and R=0.64 for the bottom) were found in the smallest, southernmost Southern California Current System (SCCS) region. The relatively complex coastline and limited resolution of island chains in this region (Fig. S2) may contribute to this decreased skill relative to other regions, but the correlation for monthly anomalies remains >0.6 even in this most challenging of systems. The SCCS bottom temperature similarly exhibits the lowest KGE (0.724), attributable to both lower correlation and variance relative to that seen in GLORYS. This was similarly the case for the bottom temperature in British Columbia (BC), which was the only other region with a KGE below 0.8.

Matching satellite-derived chlorophyll climatologies and time series (Fig. 18) proved more challenging than temperature. The monthly chlorophyll climatologies had moderate (R≥0.8 NCCS, CCCS) to high (R≥0.9, GOA, BC, SCCS) consistency with OC-CCI-based estimates for all systems except the Bering Sea (Fig. 18, left column). In the Bering Sea, NEP10k has a pronounced late-spring-to-summer peak approaching 4 mg Chl m⁻³, while OC-CCI estimates comparable intermediate concentrations of ∼2 mg Chl m⁻³ for all months but January and December. Similar, though less marked, discrepancies were found in the Gulf of Alaska. In the California Current, chlorophyll concentrations in both NEP10k and OC-CCI peak in the late spring and summer, consistent with the timing of the upwelling season. NEP10k estimates tend to drop more rapidly than OC-CCI estimates in the fall, with the central CCS exhibiting a secondary fall peak not found in NEP10k. Notably, the shelf chlorophyll comparisons in Fig. 18, which focus on temporal chlorophyll variability within a defined region, are not log transformed. This amplifies the discrepancies at the higher end of the observed range relative to those in the full domain, which focus on the model's ability to capture order-of-magnitude cross-ecosystem differences (Fig. 10).

Figure 18 Regional chlorophyll time series comparisons. Regional shelf (<500 m) surface chlorophyll monthly climatologies (left column) and anomaly time series (right column) for the satellite-derived OC-CCI product (orange) and NEP10k (green). Pearson correlation coefficients are reported for both climatologies and anomalies; background shading in the monthly anomaly time series plots indicates the oceanic Niño index produced by the NOAA Climate Prediction Center for context.

KGE for full chlorophyll time series are more moderate than those achieved for temperature, though only the GOA and the SCCS exhibit values below 0.6. For the GOA, this value was attributable to both the relative bias and variance, while in the SCCS, lower correlation and relative variance were the primary contributors to a lower KGE. Regional monthly anomaly time series for NEP10k chlorophyll were generally weakly correlated with OC-CCI (Fig. 18, right column), with most R values slightly below 0.4. While these correlations are significant (p<0.01), their modest values temper expectations for actionable chlorophyll forecasts. A possible exception is found in the northern California Current, where high correlation (R=0.58) provides some ground for optimism. Conversely, simulated and OC-CCI chlorophyll anomalies in the Bering Sea were uncorrelated ($R=-\mathrm{0.01}$). We emphasize that interpretation of both NEP10k's correspondence and misfits in Fig. 18 must be moderated by uncertainties associated with the derivation of satellite-based ocean color products in coastal waters.

3.2.1 Bering-Sea-specific indicators

As discussed in Sect. 1, the eastern Bering Sea has one of the most prolific demersal/benthic fisheries in the world, and its ecosystem dynamics are strongly shaped by fluctuating seasonal sea ice. Compared to the trawl results, NEP10k trawl-equivalent bottom temperature (Fig. 19) in the Bering Sea tends to be biased slightly warm, particularly in the mid-shelf region that approximately corresponds with the area of maximum/minimum September ice edge extent reported by Wang et al. (2014). The model exhibits a modest cold bias, in contrast, on the inner shelf of the southeastern Bering Sea. The NEP10k model, however, robustly reproduces interannual variability of the cold pool area (CPA) indices, with the best performance at the higher temperature thresholds (Fig. 20). The model does tend to underrepresent the CPA delineated by the coldest threshold (water temperature $\le -\mathrm{1}$ °C; dark blue in Fig. 20). This is emphasized by a particularly low KGE (−0.111), which is due to a particularly low relative bias and high relative variance. KGE for the $\le -\mathrm{1}$ °C threshold is also low, similarly due to relative variance and bias. However, the correlation remains high across thresholds, and there is minimal bias at the higher thresholds (i.e., water temperature ≤1 or 2 °C; lighter blues in Fig. 20). Critically, the simulation captures the very small CPAs in recent years that have been linked to recent declines in the lucrative snow crab fishery (Szuwalski et al., 2023).

Figure 19 Bering Sea cold pool extent. Comparison with AFSC Bering Sea summer trawl. Marker size is scaled by the number of annual data points that comprise the mean. The color map in a and b emphasizes the 2 °C transition point for consistency with the threshold value for identifying the cold pool. The black outline delineates the southeastern Bering Sea; trawl data collected from this region are used to calculate the Bering Sea summer cold pool extent and index.

Figure 20 Southeastern Bering Sea cold pool area index. Comparison of the cold pool index time series derived from the AFSC bottom trawl survey data (a) and the spatially and temporally consistent NEP10k bottom temperature output (b) following the methods described in Rohan et al. (2022) and the AFSC cold pool software repository. The plots report the fraction of the total survey of the southeastern Bering Sea trawl area (outlined in the figure above) that exhibits bottom temperatures under the specified thermal thresholds. We report Spearman correlation values between NEP10k and trawl indices in the bottom panel (b).

The NEP10k simulation does overestimate the sea ice concentration, particularly in the northern Bering Sea (Fig. 21). However, the contours for the 10 % and 50 % sea ice concentrations correspond with observations fairly well from January through April, suggesting that the simulation generates a reasonable spring sea ice extent. NEP10k ice extent time series for the southeastern Bering Sea (Fig. S16 in the Supplement) are highly correlated with the satellite product, though NEP10k does overestimate the coverage area, which may be consistent with the ∼0.5 °C Bering Sea cold bias noted in Fig. 2.

Figure 21 Bering seasonal sea ice concentration and spatial extent. Comparison of spatial patterns in the Bering Sea monthly mean NEP10k sea ice concentration against NASA satellite estimates (Cavalieri et al., 1996). Black contours indicate the positions of the 10 % and 50 % sea ice concentrations.

3.2.2 Gulf-of-Alaska-specific indicators

NEP10k successfully simulates the two leading localized modes of SSH variability identified by Hauri et al. (2024) that can predispose the Gulf of Alaska to extreme physical and biogeochemical events (Fig. 22). The first two principal components (PCs) of the empirical orthogonal analysis of monthly NEP10k SSH in the Gulf of Alaska have spatial patterns that are consistent with the CMEMS SSH product, with significantly correlated spatial loading patterns in both cases (EOF1 R=4.2, EOF2 R=0.95; Fig. 22, top panels). The NEP10k-generated NGAO and GOADI time series are also in good agreement with satellite altimetry observed over the corresponding region and time frame, particularly at lower frequencies (Fig. 22, bottom panels). These two modes of variability comprise 47 % and 34 % of the variance in the model and observed SSH, respectively, suggesting that they may be somewhat overprominent in the model relative to other sources of SSH variability.

Figure 22 GOA SSH EOFs and principal component time series. Spatial maps of the first (top row) and second (middle row) EOFs for satellite (left) and NEP10k (right) SSH variability. These are complemented with time series comparisons (monthly, left; 6-year running mean, right) for the first two principal components (NGAO, top row; GOADI, bottom row) from the empirical orthogonal function analyses of Gulf of Alaska sea surface height for NEP10k (orange) and the CMEMS satellite product (navy). R values indicate the Pearson correlation coefficient calculated between NEP10k and the satellite product, all of which are significant at p<0.001. X axis labels indicate 1 January of the specified year.

Composites of environmental conditions when the second PC, the GOADI, is below or above 1 demonstrate the impact of downwelling and relaxation of downwelling conditions, respectively, on shelf habitat in the Gulf of Alaska (Fig. 23). Relaxation of downwelling is associated with colder, lower-oxygen and more-acidic shelf waters from the enhanced intrusion of deep water. Conversely, positive phases of the GOADI exhibit significantly warmer bottom temperatures and elevated levels of bottom dissolved oxygen and aragonite saturation state.

Figure 23 GOADI composites. Composites of important ecological conditions during the positive (GOADI>1; 44 months out of 324) and the negative ($\text{GOADI}<-\mathrm{1}$; 45 months out of 324) phases of the Gulf of Alaska Downwelling Index (GOADI). Grid cells are colored where the composite differs significantly from 0 (student t-test, p<0.05).

3.2.3 California-Current-specific indicators

Seasonal upwelling plays an important role in CCS ecosystem dynamics, having bottom-up driving effects on primary productivity in this eastern boundary upwelling system (Sect. 1, Jacox et al., 2016). Summer upwelling conditions are evident in the map of vertical velocity (Fig. 24) with, on average, a predominantly positive/upward signal across the approximate mixed layer depth (30 m) from March through August, similar to that reported in Jacox et al. (2018). Monthly climatologies of NEP10k simulated vertical transport across 30 m demonstrate high correlation with the Jacox et al. (2018) CUTI metric, with R values above 0.92 at representative latitudes (Fig. 24). Correlations between the Jacox et al. (2018) monthly CUTI anomaly time series and corresponding NEP10k vertical transport are also significant, but the relationship is strongest at more northern latitudes (R=0.76 at 45° N) and drops off at more southerly latitudes (R=0.30 at 35° N). The 33N bin also exhibits the lowest KGE (0.248) for the full CUTI time series comparison, due to both a relatively low Pearson correlation coefficient and a fairly high relative variance. It is important to note, however, that the NEP10k and the ROMS model in Jacox et al. (2018) are forced by different atmospheric reanalysis products; thus, it may not be surprising that they differ in high-frequency variability. Additionally, the differences in methodologies such as approximating using a constant reference depth of 30 m for NEP10k could contribute to departures.

Figure 24 CCS upwelling indices. Spring/summer (March–August) vertical velocity (map) at a depth of 30 m. Bins of 1° are indicated by the black outline and used for integrating vertical transport. This (blue line) is compared against the Jacox et al. (2018) ROMS CUTI metric (orange line) at several latitudes, decomposing the time series into monthly climatology (left) and anomalies (right). Pearson correlations (R) are reported in the upper-right corner of each time series panel; all correlations are significant (p<0.001).

NEP10k trends in dissolved oxygen reproduce offshore CalCOFI trends (Fig. 25), with strongest declines occurring at around 300 m and becoming less pronounced with depth. In the California Bight, however, NEP10k exhibits positive trends (most pronounced at a depth of 100 m) where the CalCOFI time series exhibit declining trends in dissolved oxygen levels. Many of the stations exhibiting discrepancies in the NEP10k are not statistically significant (p<0.05), and it should be noted that some of the time series are quite variable, with linear trends being sensitive to the time frame analyzed. Indeed, direct point-to-point comparisons against the CalCOFI bottle sample dataset (Figs. S26–S29 in the Supplement) demonstrate that, while NEP10k broadly reproduces temperature and salinity (Fig. S26, r≥0.89) and biogeochemistry (e.g., Fig. S27, r≥0.96) conditions in the southern California Current System, the model was more challenged to represent the temporal variation observed across decades for individual sampling sites and depth strata (Fig. S28). Agreement was best at the surface and for temperature but generally decreased with depth. Skill improved when values averaged across the CalCOFI sampling grid were considered (Fig. S29).

Figure 25 CCS trends in dissolved oxygen at CalCOFI stations. Linear trends in subsurface dissolved oxygen (O₂) at CalCOFI stations for NEP10k (left) and the CalCOFI dataset (right) calculated over the time frame of the NEP10k hindcast (1993–2019). Black markers indicate where station trends are significant (p<0.05), following Bograd et al. (2008).

3.3 Computational performance and scalability

As described in Sect. 1, the goal of the NEP10k configuration is to provide a simulation capable of skillfully resolving fisheries-critical features with manageable computational cost to allow for ensemble predictions and projections. Our baseline simulation averaged just over 5.3 h of wall clock time per hindcast year while distributing the 342×816 grid (cross-shore × along-shore) across a 32×80 decomposition (Fig. 26, green circle) and using a 400 s baroclinic time step and a 1200 s thermodynamic and tracer time step. After land masking, the run used 2036 PEs, yielding roughly 10 800 PE hours per simulation year on the c5 partition of NOAA's GAEA supercomputer. The 27-year hindcast produced herein thus requires ∼292 000 PE hours, while 1200 years of retrospective seasonal forecasts (e.g., Ross et al., 2024) would require approximately 13 million PE hours.

Figure 26 Computational scalability efficiency. Amount of computer wall clock time used for completing 1 year of NEP10k simulation with a given number and configuration of processing elements (PEs). Markers indicate a given simulation's PE decomposition for diving in the horizontal model domain prior to omitting PEs that do not contain any ocean grid cells. The diagonal lines indicate constant computational cost (processes × time) relative to the 40×40 (blue square) reference simulation. The two hollow markers represent simulations wherein the thermodynamics time step was set to a dynamics time step (i.e., reduced from 1200 to 400).

The NEP10k computational cost is comparable to the recently published Northwest Atlantic regional MOM6 configuration (NWA12) of Ross et al. (2023), which used a 40×40 layout (1200 PEs after land masking) to generate 1 simulation year in about 9 h (about 10 800 PE hours per simulation year). While NWA12 was a larger domain, NEP10k required smaller baroclinic and thermodynamic time steps for stability (400 versus 600 s and 1200 versus 1800 s, respectively). The instability at longer time steps in the NEP10k configuration primarily occurred in the vicinity of the Aleutian island chain, where strong currents could be generated within tight channels.

Computational scaling tests showed that increases in throughput were achievable, but returns fell considerably below the ideal 1:1 scaling between the processor count and the wall clock time (Fig. 26). An approximate doubling of PEs from 2038 to nearly 4000, for example, only decreased the wall clock time for a simulation year from ∼5.3 to ∼4.2 h (compare the green circle and the purple diamond in Fig. 26). The decreased scaling is not unexpected, as higher processor counts decompose the model grid into increasingly granular tiles, taxing communication across PEs. This effect can also be seen when comparing the performance of the 32×80 baseline setting, which maximizes the number of interior to exterior cells on a PE by decomposing the 342×816 grid into squares, versus the approximately 10 % slower 50×50 decomposition that relies on rectangular elements. Scaling from the base configuration to lower processor counts, in contrast, is relatively strong, supporting the viability of running simulations on smaller supercomputing systems.

Consistent with the findings of Ross et al. (2023), we found considerable computational benefit from leveraging MOM6's capacity to have a longer thermodynamic and tracer time step than the baroclinic time step (closed versus open symbols in Fig. 26). Throughput was nearly doubled when the thermodynamics and tracer time step was 3 times longer than the baroclinic time step.

4 Discussion

There were three primary design criteria for the NEP10k model. The first was that a coastwide configuration was needed to address coastwide challenges arising from climate change, such as shifting fisheries distributions across state and international boundaries. The second was that the model must resolve and accurately reproduce enough of the physical and biogeochemical drivers of ocean change in and across the disparate ecosystems within the domain to support ecosystem and fisheries applications. The third was that the model must be suited, both computationally and in terms of model skill, for ensemble predictions and projections. The comprehensive model evaluation herein suggests that the NEP10k configuration meets these design criteria sufficiently to provide a basis for initial applications and a robust foundation for further model improvement. The comparison against large-scale physical and biogeochemical patterns in Sect. 3.1 showed that a single physical–biogeochemical modeling framework could robustly capture the primary physical and biogeochemical contrasts between the EBS, GOA, and CCE (Figs. 2–5, 7–8, 12–15). Simulation fidelity extended to seasonal patterns in most quantities (Figs. 17 and 18, Figs. S1–S15) and robust matches to interannual variations for many, even within limited regions of the domain (Figs. 16, 19, 21, 24, and 25). While biases were present, and at times prominent, the skill achieved supports NEP10k's current utility. This discussion will focus on model characteristics contributing to successes and on further model developments that may ameliorate current limitations.

A central challenge for NEP10k was the representation of physical and biogeochemical processes governing a large range of ecosystems, from subtropical to polar and oligotrophic to eutrophic. Success in this regard requires model formulations and parameterizations that are robust across regimes. For ocean physics, one advance that led to notable improvement was the replacement of the submesoscale restratification parameterization of Fox-Kemper et al. (2011) with that of Bodner et al. (2023). The Fox-Kemper parameterization requires a single choice for the submesoscale front length, while Bodner diagnoses the front length from the ocean state, revealing considerable variability with season and latitude. Smaller front lengths at high latitudes proved critical to limiting deep mixing biases in the western Bering Sea, while longer front lengths further south were critical in limiting shallow mixed layer biases in the Gulf of Alaska and California Current (Fig. 4). Though the more dynamic Bodner scheme did not eliminate MLD biases, we did find that it improved them considerably relative to the Fox-Kemper et al. (2011) parameterization, where a single characteristic submesoscale frontal length scale forced one to exacerbate one bias or the other (Fig. S18 in the Supplement).

For biogeochemistry, starting with a model designed for global applications provided a sound starting point for achieving cross-system skill. Evaluation of the shelf-scale fidelity of global models, however, is generally limited by their often coarse resolution (e.g., Stock et al., 2014, 2020). A key addition to extend skill in NEP10k to coastal regions was an additional phytoplankton size class, which allowed the model to better resolve the coastal diatoms responsible for high chlorophyll concentrations along the coast. This expanded formulation was initially developed by Van Oostende et al. (2018) for use in the California Current, where it was shown to improve the resolution of both very high coastal chlorophyll concentrations and the biogeochemical signals that can be associated with them (e.g., coastal hypoxia). These benefits can be seen in the generally high coastal (relative to open ocean) chlorophyll levels along the U.S. West Coast (Fig. 10) and the robust depiction of the hypoxic boundary layer depth (Fig. 13). The most glaring chlorophyll bias is the model's tendency to underestimate winter/fall OC-CCI-estimated chlorophyll in the nearshore EBS (Fig. 10), which degrades the seasonal chlorophyll fidelity for this region (Fig. 18). Satellite-based estimates in shallow regions of the EBS actually peak during these months despite cold, dark, and vigorously mixed conditions, suggesting potential contamination of chlorophyll estimates in turbid coastal waters (Dierssen, 2010; Schofield et al., 2013). A recent study in the Arctic, for example, suggests that global satellite chlorophyll algorithms may overestimate chlorophyll by over a factor of 2 (Li et al., 2024).

Other chlorophyll and plankton misfits require additional scrutiny. The tendency to overestimate offshore spring and summer chlorophyll along the margin separating the Gulf of Alaska and the California Current, for example, may reflect biases in dust delivery, dust solubility, or iron scavenging in this iron-limited region. The relatively persistent and strong iron limitation in the offshore waters of the California Current in NEP10k, however, may already exceed the “mosaic” of alternating N and Fe limitation suggested by some prior studies (Messié and Chavez, 2015; Moore et al., 2013; Till et al., 2019). A spatially indiscriminate iron tuning is thus unlikely to resolve these biases. They may also arise, however, from misrepresented grazing controls. NEP10k's skill in simulating mesozooplankton biomass is limited to capturing first-order cross-ecosystem and seasonal biomass contrasts (Fig. 11), with the patchiness in mesozooplankton biomass in net tow data being underrepresented. There are also some systematic biases, such as the tendency for mesozooplankton populations to be displaced offshore and biased low relative to observations during the summer upwelling season in the California Current. Previous work (e.g., Batchelder et al., 2002) has suggested that zooplankton may enlist diurnal vertical migration to avoid being swept offshore, alternating between surface feeding in offshore currents at night and predator avoidance in inshore flowing currents during the day. Such behavior is not included in NEP10k but could increase mesozooplankton biomass and shift the distribution inshore.

Capturing mean spatial and seasonal patterns is a critical starting point for any model intended for ecosystem/fisheries science and management applications. Many applications, however, require the capacity to anticipate change across seasonal to multidecadal management time horizons (Tommasi et al., 2017). The robust representation of surface and bottom temperature variability (Fig. 16) provides a promising start in this regard. Temperature anomalies are a first-order indicator of ecosystem conditions and a primary determinant of habitat viability (e.g., Deutsch et al., 2015), and temperature extremes are a primary source of ecosystem stress in a changing climate (e.g., Frölicher et al., 2018). The robust representation of surface and bottom water anomalies at a regional scale and for shallower waters (<500 m), combined with the growing capacity of global prediction systems to anticipate fluctuations in large-scale climate drivers (e.g., ENSO), supports the potential viability of predictive applications. Retrospective forecast experiments are underway to assess this. NEP10k was less successful, however, in capturing coastal chlorophyll anomalies (Fig. 17). The correlation with monthly chlorophyll anomalies was only marginally significant in most systems, approaching useful levels (i.e., R∼0.6) in the NCCS. This weaker correlation was not necessarily surprising, given the volatile and patchy nature of coastal chlorophyll and observing challenges in such environments, but points to the need for further scrutiny of both the model and observations before predictive chlorophyll applications can be realized in most systems.

Possibly the most critical metrics for ecosystems and fisheries applications considered herein were the region-specific quantities considered in Figs. 18–25. These were drawn from existing management-linked documents, such as the “State of the Ecosystem” reports created by NOAA's National Marine Fisheries Service to strategically inform management decisions. Evaluations against the admittedly limited set of region-specific fisheries metrics herein was generally positive. Perhaps the most striking of these successes is the fidelity with which NEP10k reproduces the Bering Sea cold pool relative to over 2 decades of Alaska Fisheries Science Center bottom trawl data (Figs. 19 and 20). The model's representation of these metrics was improved during the course of development when an excess of shear-driven mixing on the Bering shelf was identified and addressed with an adjustment of Jackson et al. (2008) shear mixing parameterization. The addition of a simple scaling factor for the geometric limitation imposed by this formulation was found to be the most effective way to pragmatically calibrate the shear-driven mixing to better produce observed values for both mixing and bottom temperature. A more comprehensive analysis of this parameterization and its impact on Bering Sea dynamics is currently underway (Seelanki et al., 2025) and will inform regional MOM6 shear mixing parameterization for mixed turbulence regimes.

While NEP10k's overall representation of variations in the Bering Sea cold pool extent was excellent, the model did underestimate the summer extent of the coldest bottom water ($<-\mathrm{1}$ °C; darkest blue in Fig. 19). This seemingly conflicts with NEP10k's overrepresentation of seasonal sea ice extent (Figs. 20 and S16) because greater sea ice extent and coverage tend to be associated with a more extensive cold pool (e.g., Wyllie-Echeverria and Wooster, 1998). The model does achieve substantial winter levels of cold bottom water (Fig. S19 in the Supplement), but they erode more quickly than observed in May and June, just prior to the trawl season. This decline in bottom coverage by the coldest water mass category coincides with a dramatic monthly reduction in NEP10k's SEBS sea ice extent relative to satellite estimates (Fig. S20 in the Supplement, May–April and June–May). The drivers of this bias will be explored. We emphasize, however, that simulated Bering Sea ice variations in NEP10k are highly correlated with observations (Fig. S16), suggesting the potential for predictive applications despite the mean sea ice bias.

NEP10k reproduction of localized modes of low-frequency climate variability in the Gulf of Alaska (NGAO and GOADI, correlation with satellite-derived PCs >0.65, Fig. 22) holds promise for potential for multiyear to decadal fisheries applications in the GOA. These modes of variability map onto important ecosystem drivers such as bottom temperature and aragonite saturation state (Figs. 22 and S17 in the Supplement) and can contribute to extreme compound events that can have severe consequences for marine ecosystems (Hauri et al., 2024). Understanding of the relationships between SSH variability and shelf ecosystem conditions will be aided by the growing availability of physical and biogeochemical observations of GOA bottom conditions. Increasing the horizontal resolution of the NEP10k configuration may further improve representations of important regional GOA ecosystem features. For example, sea surface heights south of the Aleutian island chain, central to the Alaska Gyre, are lower than observed in reference datasets (Fig. 5) and could improve with better resolution of opposing horizontal flows, specifically the southwestward Alaska stream and eastward Subarctic or Aleutian Current. Higher resolution may also improve the representation of transports through the Aleutian island chain, which can significantly impact water mass properties in the Bering Sea (Stabeno et al., 1999).

Finally, in the California Current System, our regional assessment focused on ecosystem-critical seasonal upwelling and source water trends. NEP10k's climatological vertical transport at 30 m along the continental U.S. West Coast is highly correlated (i.e., R values≥0.93, Fig. 24) with the CUTI metric published by Jacox et al. (2018). Similarly, reproduction of multidecadal trends in dissolved O₂ (Fig. 25) observed in the CalCOFI record was an important benchmark, indicative of the model's ability to capture processes driving ecologically consequential deoxygenation in the southern CCE (Bograd et al., 2008). While these findings further support the suitability of the current NEP10k configuration for ecological applications, continued model development will seek to understand and improve localized performance. For example, warm/cold-biased climatological surface/bottom temperatures in both CCCS and SCCS (Fig. 17), underrepresentation of climatological upwelling and low correlation in upwelling monthly anomalies (33N in Fig. 24), and underrepresentation of deoxygenation trends in the Southern California Bight (200 m and 300 m depths in Fig. 25) suggest that we may not be adequately representing the physical processes that influence these conditions due to excessive stratification in the southern CCE. Given the complex bathymetry and circulation that impact these processes in the Southern California Bight (e.g., Hickey, 1992), this is another instance where increased spatial resolution may improve model performance. Indeed, at a 10 km resolution, the SCCS shelf (where depth ≤500 m) extends as little as a single grid cell (or less) offshore. However, while higher-resolution (i.e., ∼5 km) simulations are currently underway, any benefits of doubling the horizontal resolution will need to be balanced against the roughly 8-fold increase in computational cost (i.e., 2-fold for each horizontal dimension and an additional 2-fold increase for the need to shorten the time step needed to maintain Courant–Friedrichs–Lewy stability). For applications wherein many-fold higher resolution is necessary, it may be more practical to utilize a smaller, higher-resolution nested domain (e.g., modeling the Salish Sea in Khangaonkar et al., 2018) that can be forced by the NEP10k at the open boundaries, rather than increasing the resolution for the full NEP10k domain. Continued NEP10k development will incorporate comparison against a broader array of local observation datasets similar to that of CalCOFI. Such extensive observation records are invaluable for better understanding and evaluating model performance, particularly in regions that may not be well represented in relatively coarse, gridded data products. However, it is important to approach such comparisons with realistic expectations. As shown in Fig. S28, NEP10k poorly reproduces temporal variability (i.e., low Pearson correlation coefficients) of repeated samplings of individual stations across multiple years. This is not surprising because the NEP10k hindcast does not assimilate observations, and thus any biases in the mean locations of fronts and other features are compounded by stochastic mesoscale and submesoscale features whose precise locations and timing will not match those observed. Indeed, more coherent patterns emerge after averaging over such features (e.g., Figs. S29, 17, and 20), which demonstrates that NEP10k's strength and utility are in representing reasonable approximations of ecologically important environmental conditions rather than an exact reproduction of in situ observations.

5 Conclusions

The results presented herein demonstrate that NEP10k is “fit for purpose” – in terms of both model skill and computation cost – for numerous living marine resource management applications across multiple time horizons. The model also establishes a basis for community evaluation to assess against a much broader set of fisheries and ecosystem metrics and a basis for co-development with fisheries scientists and managers to address identified limitations and maximize model utility. As part of NOAA's Changing Ecosystems and Fisheries Initiative, the community contributing to this effort has grown tremendously, facilitated by the open development of MOM6, COBALT, and preprocessing and analytical scripts made available via the CEFI GitHub. With increasing input from collaborators and co-development with end users, ongoing model development will prioritize NEP10k representation of key ecosystem indicators to maximize the utility of climate change projections and forecasts for living marine resource management.

Appendix A

Table A1 Notable parameters, their current names, and the associated values used in the physical ocean (MOM6) component of the model and relevant references. BGC denotes biogeochemistry; SAL denotes self-attraction and loading. Bold text indicates where parameter choices differ from Ross et al. (2023). Comprehensive documentation of physical MOM6 parameters can be found in MOM_parameter_doc.all (Drenkard et al., 2024a).

Table A2 Ocean diagnostics used for evaluating the NEP10k hindcast.