Species distribution models (SDMs) have become key tools for describing and predicting species habitats. In the marine domain, environmental data used in modeling species distributions are often remotely sensed, and as such have limited capacity for interpreting the vertical structure of the water column, or are sampled in situ, offering minimal spatial and temporal coverage. Advances in ocean models have improved our capacity to explore subsurface ocean features, yet there has been limited integration of such features in SDMs. Using output from a data-assimilative configuration of the Regional Ocean Modeling System, we examine the effect of including dynamic subsurface variables in SDMs to describe the habitats of four pelagic predators in the California Current System (swordfish Xiphias gladius, blue sharks Prionace glauca, common thresher sharks Alopias vulpinus, and shortfin mako sharks lsurus oxyrinchus). Species data were obtained from the California Drift Gillnet observer program (1997-2017). We used boosted regression trees to explore the incremental improvement enabled by dynamic subsurface variables that quantify the structure and stability of the water column: isothermal layer depth and bulk buoyancy frequency. The inclusion of these dynamic subsurface variables significantly improved model explanatory power for most species. Model predictive performance also significantly improved, but only for species that had strong affiliations with dynamic variables (swordfish and shortfin mako sharks) rather than static variables (blue sharks and common thresher sharks). Geospatial predictions for all species showed the integration of isothermal layer depth and bulk buoyancy frequency contributed value at the mesoscale level (\u3c 100 km) and varied spatially throughout the study domain. These results highlight the utility of including dynamic subsurface variables in SDM development and support the continuing ecological use of biophysical output from ocean circulation models
