Dynamic occupancy modeling of temperate marine fish in area‐based closures

Abstract Species distribution models (SDMs) are commonly used to model the spatial structure of species in the marine environment, however, most fail to account for detectability of the target species. This can result in underestimates of occupancy, where nondetection is conflated with absence. The site occupancy model (SOM) overcomes this failure by treating occupancy as a latent variable of the model and incorporates a detection submodel to account for variability in detection rates. These have rarely been applied in the context of marine fish and never for the multiseason dynamic occupancy model (DOM). In this study, a DOM is developed for a designated species of concern, cusk (Brosme brosme), over a four‐season period. Making novel use of a high‐resolution 3‐dimensional hydrodynamic model, detectability of cusk is considered as a function of current speed and algae cover. Algal cover on the seabed is measured from video surveys to divide the study area into two distinct regions: those with canopy forming species of algae and those without (henceforth bottom types). Modeled estimates of the proportion of sites occupied in each season are 0.88, 0.45, 0.74, and 0.83. These are significantly greater than the proportion of occupied sites measured from underwater video observations which are 0.57, 0.28, 0.43, and 0.57. Individual fish are detected more frequently with increasing current speed in areas lacking canopy and less frequently with increasing current speed in areas with canopy. The results indicate that, where possible, SDM studies for all marine species should take account of detectability to avoid underestimating the proportion of sites occupied at a given study area. Sampling closed areas or areas of conservation often requires the use of nonphysical, low impact sampling methods like camera surveys. These methods inherently result in detection probabilities less than one, an issue compounded by time‐varying features of the environment that are rarely accounted for marine studies. This work highlights the use of modeled hydrodynamics as a tool to correct some of this imbalance.


| INTRODUC TI ON
The use of species distribution modeling as a tool for scientists and environmental managers has seen a substantial increase in the last three decades, driven by a growing demand for knowledge of species' ranges and facilitated by increasingly powerful computing resources (Barbosa & Schneck, 2015). Species distribution models (SDMs) mathematically represent the relationship between species records and features of the environment, often with the intention of predicting suitable ranges for the target species (Franklin, 2010).
Depending on the aims of an investigation and the type of data collected, several approaches are available to researchers. For many applications where surveys have been planned in advance of statistical analysis, standard methods such as generalized linear models, generalized additive models, and random forest are frequently used to model occurrence and abundance data (Elith et al., 2011).
However, given the time and cost of performing a systematic survey, especially in the marine environment, researchers often use archived datasets to perform modeling investigations (Araújo & Guisan, 2006). In many cases, these datasets only contain records on species presence, resulting in a situation where no absence data are available (Elith et al., 2011). For this reason, presence-only models such as MaxEnt (Phillips, Anderson, & Schapire, 2006) have been developed that do not require absence data, but rather generate a large number of pseudo-absences from the study area (Phillips et al., 2006).
While each approach outlined above has advantages, both fail to address the issue of detectability (Monk, 2014). Both approaches assume that detection probability is invariant, that is, that the target species is perfectly observed whenever it is present (Yackulic et al., 2013). This is often not the case when studying cryptic species and is an important consideration in the marine environment when sampling methods often do not result in direct observation of the study environment, for example, when using trawl or camera surveys (Monk, 2014). In a review of 108 articles that used MaxEnt, Yackulic et al. (2013) found that only 14% mentioned detection probability.
This failure to address detectability introduces error to estimates of occurrence for the species being modeled and can result in erroneous reporting of covariate effects (Guillera-Arroita, Lahoz-Monfort, MacKenzie, Wintle, & McCarthy, 2014).
Site occupancy modeling (SOM) allows occupancy and detectability to be analyzed hierarchically as two separate processes, accounting for the problem of imperfect detection (MacKenzie et al., 2002). Monk (2014) provides a good overview of the need to adopt this class of model in the marine environment; while available for a similar period as MaxEnt, the SOM has received much less attention (e.g., Coggins, Bacheler, & Gwinn, 2014) in marine ecology investigations. The multiseason occupancy model (or dynamic occupancy model for its Bayesian counterpart (DOM)) allows occupancy, detection, local colonization, and extinction to be accounted for across several sampling seasons (MacKenzie, Nichols, Hines, Knutson, & Franklin, 2003). Seasons refer to primary sampling periods within which the population is assumed closed but between which the population can be subject to local extinction and colonization. The model assumes that at least one site has been visited more than once within a sampling period and that the true occupancy state is imperfectly observed, that is, it is a latent variable of the model. As such, the model can be thought of as a nonstandard GLMM with a binary random effect equal to 1 where the site is occupied by the target species and 0 where it is not (Kéry, 2010). Each of the four probabilities within the DOM (initial occupancy, colonization, extinction, and detection) can be modeled as a function of a set of covariate data or set as constant across sites within a given sample period (MacKenzie et al., 2003;Royle & Kéry, 2007). Where covariates are used to model probabilities, these must represent variability in the environment along temporal scales relevant to the phenomenon under study. Extinction and colonization effects, for example, require seasonally varying covariates. By contrast, detection effects require covariates that vary over much shorter temporal scales, allowing differences in detectability to be discerned within relatively short sampling periods. One possible reason why these models have not received as much attention in the marine environment as they have for terrestrial studies is the prohibitive cost of marine sampling.
In order to satisfy the assumptions of the dynamic occupancy model, repeat visits are required within each of a number of seasons to collect both response and environmental data.
Covariates used for modeling that have been collected in situ at the time of making species observations are often assumed to better describe observed patterns in species distributions (Franklin, 2010).
Recent research has shown, however, that this is not always the case and that modeling studies can benefit from a combination of in situ sampling and remote sensing data (Niedballa, Sollmann, Mohamed, Bender, & Wilting, 2015). Moreover, Newton-Cross, White, and Harris (2007) demonstrated that remotely sensed or computer generated data can be more effective than data collected in situ for accurately predicting the occurrence of some terrestrial based species.
Again, this is important in the marine environment where sampling is often expensive and time-consuming compared to terrestrial studies. As such, marine researchers often have to rely on remote sensing (Brown, Sameoto, & Smith, 2012) and modeled data (Rattray, Ierodiaconou, & Womersley, 2015) to supplement in situ sampling.
These latter datasets often come in the form of hydrodynamic models that mathematically represent tidal and wave forcing of the marine environment (Gunn & Stock-Williams, 2013). While widely used in engineering and physical oceanography (e.g., Chen et al., 2011;McMillan & Lickley, 2008), their use in marine ecological investigations has been limited. This is likely due to the expense, in terms of time and computational power, of setting up a hydrodynamic model that accurately reflects conditions at spatial and temporal scales relevant to ecological processes and organism behavior.
This investigation aims, for the first time, to create a dynamic occupancy model to demonstrate the effectiveness of such an approach for a temperate marine fish species in a closed fishing area.
The study generates unbiased estimates of occupancy and compares these to occupancy estimates obtained by survey methods alone. Inputs to the dynamic occupancy model comprise observed, derived, and simulated data. Observed data include observations of the target species, algae cover, and geomorphological complexity from video surveys; derived data include depth measured using multibeam echosounder (MBES) data, and terrain attributes derived from the MBES depth data. Additionally, the study makes novel use of current velocities simulated using a high-resolution hydrodynamic model to demonstrate the utility of including these data in fisheries monitoring investigations.

| Analysis overview
This study takes a multifaceted approach to producing a dynamic occupancy model (DOM) for a species of temperate marine fish at a remote rocky outcrop in the central Gulf of Maine. The model was specified in a Bayesian framework to account for a small amount of separation in the detection data and to allow the use of a finite sample estimator that generates more accurate estimates of occupancy for a small sample size (Royle & Kéry, 2007). The analysis used data observed or derived from two primary datasets: a series of noninvasive underwater video surveys collected over four sampling seasons and a MBES survey conducted for the Gulf of Maine mapping initiative (SAIC, 2005). Observations of the response, cusk, were made using the video footage. During each video survey, detections of cusk were recorded along with the time of the observation.
Simultaneous observations of algae cover and morphological complexity were made from the video footage, to be later used as explanatory variables. Separately, the digital elevation model extracted from the MBES data was used for three purposes: (a) to derive terrain attributes to assess morphological complexity over the study site; (b) to create a surface of algae cover over the study site based on empirical extinction depths; and (c) as an input to a standalone hydrodynamic model. The hydrodynamic model was used to estimate bottom current conditions at the study site, from which two outputs were generated (a) point estimates of time-varying bottom current speed at each of the video locations surveyed for cusk and (b) time-varying surfaces of current speed for the entire study area.
Once all primary data had been processed, the DOM was created.
Model coefficients for the significant terms in the detection submodel of the DOM (observed algae cover and bottom current point estimates) were obtained and used to generate the final outputs.
This was achieved by predicting the model coefficients over the algae cover surface and bottom current surfaces created above, creating spatio-temporally varying predictions of detection probability for the study area. Each step of the analysis is described in full in the following sections ( Figure 1).

| Candidate species and area
The candidate species for this investigation is cusk (Brosme brosme, Lotidae); a cryptic, bottom dwelling species found in the eastern and western Atlantic Ocean and designated a species of concern by the F I G U R E 1 Flow diagram for analyses described in the methods. Solid lines represent the flow of information to create the dynamic occupancy model, and dashed lines represent the flow of information for producing the output probability surfaces in Figure 6 National Oceanic and Atmospheric Administration (NOAA) National Marine Fisheries Service (NOAA 2009). In the western Atlantic, cusk are found from Nova Scotia in Canada to New Jersey in the USA and typically stay in deeper waters (>100 m) in these areas. Within the Gulf of Maine, however, they are typically found in shallower water owing to the relatively shallow depths of the internal Gulf (Bigelow & Schroeder, 1953). While relatively little is known about their specific life history and ecology (Davies & Jonsen, 2011), it has been noted that cusk prefer structured habitat and use kelp forests, boulder piles, and rock crevices as refugia (Auster & Lindholm, 2005;Hare et al., 2012). In addition, they are considered to be weak swimmers (Bigelow & Schroeder, 1953), so make an ideal species for study on how their behavior is affected by variability in movements of the water column. The species also has a small home range (Dultz, 2013), making it suited to the assumptions of the dynamic occupancy model.  (Sherwood & Grabowski, 2016) and supports large resident populations of several commercially important fish species (Grabowski, McGonigle, & Brown, 2010). The site displays a tripartite zonation of macroalgae around the summit with each of the three zones reaching record depths for boreal-subarctic waters: leathery macrophytes to 40 m, foliose red algae to 50 m, and crustose algae to 63 m (Vadas & Steneck, 1988). The Ledge comprises morphologically complex granite shoaling at 10 m water depth, with sand and gravel deposits appearing around 60 m and silt dominated habitats below 80 m. To the east and west, the site is flanked by relatively deep basins (<220 m) dominated by sands, fine silt, and clay (Uchupi & Bolmer, 2008). During each sampling season, a maximum of three replicate surveys were conducted at each site, with an average of two. Camera units were deployed and left in situ for up to 1.5 hr during each survey, recording the time in and position of the camera. The camera was mounted on the sampling equipment such that the field of view was parallel to the seabed; no directional controls were in place, so the azimuthal direction varied between samples. For full camera set up see Grabowski et al. (2010). Any samples where the camera equipment landed with the camera facing into the water column and therefore unable to view the seabed were discarded and not used in any further analyses. After checking each sample for positional accuracy, videos were examined noting the time on the video any cusk were observed. For videos where cusk were present, the video time was combined with the survey start time to obtain the exact time of the observation. Where cusk were absent in a video, a time was randomly sampled from the length of the video to obtain a time for the null observation.

| Data collection and covariate generation
In addition to the observations of cusk, the videos were used to qualitatively assess fine-scale morphological complexity (high, moderate, low) and algae density. Algae density was used to classify the study area into areas of two different bottom types according to algae cover, henceforth bottom type: areas with canopy forming species and areas with no canopy forming species. Sampling effort was defined simply as the length of bottom time in each video.
Depth for each observation was obtained from 5 m resolution MBES data collected for the Gulf of Maine Mapping Initiative (SAIC 2005).
Morphological complexity was derived from the MBES data using the relative deviation from the mean value (RDMV) as recommended by Lecours, Devillers, Simms, Lucieer, and Brown (2017) in a 3 × 3 cell moving window. A breakdown of the mean, minimum, and maximum values for each of the depth strata sampled is in Table 1.

| Hydrodynamic model
Hydrodynamics were assessed using a coupled wave-current model produced using MIKE by DHI (Danish Hydraulic Institute, 2014).

The current model solves the three-dimensional incompressible
Reynolds averaged Navier-Stokes equations, while the wave model solves a fully spectral wind-swell formulation. The models are coupled to include wave-current interactions and are solved using a finite volume method over a flexible mesh that allows higher resolution in areas of interest (Danish Hydraulic Institute, 2014).
The domain for the model (Figure 2b) incorporated the Gulf of Maine from Tor Bay to Rhode Island and extends seaward off the continental shelf to allow the Gulf to respond freely to tidal forcing (McMillan & Lickley, 2008). Forcing was supplied to the model as spatially and temporally varying surface elevation from the DTU10 0.125° global tidal model (Cheng & Andersen, 2010). Calibration of the model was achieved by adjusting the value of bed resistance over a series of 13-month simulations (one-month warm-up, 12month usable data). After each calibration simulation, harmonic analysis was conducted for 67 sites within the Gulf of Maine for comparison against empirical data from Moody et al. (1984). Where disagreement between known and modeled data was unacceptably large, values of bed resistance were iteratively adjusted to fine-tune the harmonics.
Once the model harmonics were calibrated to maximally correspond to empirical data, the computational mesh was refined to increase resolution at Cashes Ledge. Here, the maximum horizontal resolution was 135 m with an average horizontal resolution of 220 m.
In addition to refining the model mesh, atmospheric forcing was introduced to the model using the National Centres for Environmental Prediction (NCEP) Climate Forecast System Reanalysis (CSFR) 6-hourly, 0.5° global weather model (Saha et al., 2010(Saha et al., , 2011. Model validation was subsequently conducted on the refined mesh model by hindcasting periods of time not included in the calibration models. Validation included the assessment of model outputs against measured wind, wave, and current data. The model was considered validated when outputs were minimally different to measured data.
All model runs were performed with a time step of 20 min to allow validation against measured data.
Full details of the model setup, calibration, and validation can be found in the Supporting Information Appendix S1. Once the model had been validated, current speeds were extracted from the waterseabed interface layer to capture current variations on the bottom.
Times for the observations of cusk were then matched to the temporally closest value of current speed. Visualization of temporal variability in current magnitude and direction was achieved using a tidal ellipse created using MATLAB (Xu, 2002). The ellipse was derived for a point in open water 100 m from the summit in order to capture the movement of water without influence from local topography.

| Dynamic occupancy model
The model was specified using notation from Royle and Kéry (2007) where ψ t is the probability of initial occupancy at time period t, ϕ t is the probability that a site remains occupied between t and t + 1, γ t is the probability that a site is colonized between t and t + 1, and p is the probability of detection. During data exploration boxplots revealed a high degree of collinearity between RDMV and fine-scale complexity observed in the video data. Fine-scale complexity was excluded from further analyses; it is a categorical variable so including it would require estimation of more model parameters. Covariates tested for ψ t were depth, fine-scale morphological complexity, RDMV, maximum current speed throughout the sampling season (CurMax) and CurMax 2 . Depth, fine-scale complexity, and RDMV as several studies have shown the importance of depth and seabed roughness for cusk habitat selection (Davies & Jonsen, 2011;Hare et al., 2012;Knutsen et al., 2009); CurMax and CurMax 2 to test for any limiting effect of hydrodynamic forcing on cusk habitat selection. Covariate effects were not included for ϕ t or γ t and were therefore assumed to be the same for all sites within each sampling period (MacKenzie et al., 2003;Royle & Kéry, 2007).
Covariates included for p in the detection submodel were current speed, sampling effort, and bottom type. Current speed as cusk are weak swimmers whose movement is likely to be affected by movements of the water column; sampling effort as the longer the camera is in the water, the more likely it is that a fish will be observed in any given sample; bottom type as canopy forming species are likely to obscure vision and therefore affect detectability of cusk. Covariate effects for each submodel were transformed using a logit link function, and all continuous covariates were standardized before analysis. In order to simplify model specification two models were assessed, one without (A) and one with (B) covariates for ψ, both of which included covariates for p.
Specification for ψ in model B: where β i are the regression coefficients, RDMV is the relative deviation from the mean value, and current max is the maximum current speed at the site in any of the sampling seasons.

Specification for p in models A and B:
where current is the current speed at the time of sampling for cusk, none is the bottom type with no canopy forming algae, and Current|none is the interaction between current speed and the bottom type with no canopy forming algae.

| Model validation and outputs
Goodness of fit (GOF) for each model was assessed using Bayesian ). The MCMC algorithm was set up with 5 chains, each sampling 50,000 draws. The first 10,000 draws were discarded as a burn-in period, and every tenth sample was stored for analysis. Chains were assessed visually for mixing and autocorrelation, and convergence was assessed using the Gelman-Ruben diagnostic (Gelman & Rubin, 1992) with values less than 1.1 considered converged. Occupancy was determined using a finite sample estimator that is recommended for a small number of nonrandomly selected sites (Royle & Kéry, 2007). This estimator allows a distinction to be made between estimates derived for parameters of the population and those derived for sites in the actual sample. The finite sample estimator reduces the variance of the point estimates generated for the sampled sites and is an additional benefit of fitting the model in a Bayesian framework (Royle & Kéry, 2007).
Continuous surfaces of modeled current speed were generated from the validated hydrodynamic model to visualize spatio-temporal variability of current speed along a tidal cycle at Cashes Ledge.
Predictions of bottom type were made based on the extinction depths of various algae previously reported at the Ledge (Vadas & Steneck, 1988). Surfaces of detection probability were then generated for the entire Ledge for the tidal cycle by predicting the occupancy model coefficients over the current speed and bottom type surfaces.
All statistical analyses and predictive mapping were carried out in the R environment (R Core Team 2015), and Bayesian inference was conducted in JAGS (Plummer, 2009).

| Cusk observations and data exploration
Of 112 replicate surveys analyzed, 39 contained cusk and 73 did not.

| Model validation and outputs
Bayesian p-values were 0.58 for model A (  (Table 2).

| D ISCUSS I ON
The problem of bias introduced by failing to account for detectability of mobile fish species when estimating their occupancy was ad-

| Cusk behavior
Cusk are slow moving weak swimmers (Bigelow & Schroeder, 1953;COSEWIC 2003), and it is therefore not unreasonable to expect them to be influenced by hydrodynamic conditions. Outputs of model A show that, indeed, cusk detectability is affected by changes in current speed along a tidal cycle. Within the noncanopy forming regions of Cashes Ledge, the increase in detectability with increasing current speed can be thought of as a proxy for increased activity of cusk. This increase in activity is interpreted in one of two ways, either as cusk searching for morphologically complex environments as refuge, or as cusk using the movement of water as an opportunity to forage for food. Cusk have been observed to prey on cunner (Tautogolabrus adspersus; Auster & Lindholm, 2005), and cunner in turn have been observed to forage more on exposed surfaces near refugia with increasing current velocity (Auster, 1988(Auster, , 1989. Other than these prey, little is known about cusk diet in the western Atlantic. In European waters, however, stomach contents analysis shows that cusk will eat a range of crustaceans and molluscs, both of which have been found in abundance in the shallower kelp dominated regions at Cashes Ledge (Vadas & Steneck, 1988;Witman & Sebens, 1992). Within this canopy forming algae region, as current speed begins to increase, detection probability decreases quite rapidly. While marine fish have been shown to use kelp habitats as both refuge and foraging grounds (Holbrook, Carr, Schmitt, & Coyer, 1990;Uhl, Bartsch, & Oppelt, 2016), the presence of canopy forming species of algae will encumber detectability, especially in relatively high flows or with high energy wave conditions (Rattray et al., 2015).
TA B L E 2 Parameters estimated for coefficients for occupancy model A (without covariates for ψ) and model B (with covariates for ψ) Note. GOF is the Bayesian p-value, where values closer to 0.5 indicate better fit. Gelman-Ruben MV is the multivariate Gelman-Ruben convergence statistic, where values close to 1 mean the model has successfully converged. Current is the current speed at the time of observation of cusk, effort is sampling effort, none is regions with no canopy forming species of algae, current|none is the interaction between current and none, complexity is RDMV, current max is the maximum current speed during the sampling season.
In previous studies of cusk habitat usage, the species has been recorded at much greater depths than those observed in this study (Davies & Jonsen, 2011;Hare et al., 2012;Nye, Link, Hare, & Overholtz, 2009). One explanation of the depths observed in this study is that the fish in this region are year round residents that have become accustomed to living at comparatively shallower depths. Bigelow and Schroeder (1953) note that cusk do not often move from bank to bank, and no seasonal spawning migrations have been noted for the species (COSEWIC 2003). The spawning season for cusk in the Gulf of Maine extends from April to July (Bigelow & Schroeder, 1953), and this might explain the large growth rate in the number of sites occupied between Spring 2007 and Summer 2007 (sampling seasons 2 and 3 respectively) as fish become more active in search of a mating partner. These two time periods together take in portions of the spawning season for cusk; it is therefore possible that the later it is in the season, the more active the fish are in their search.
An assumption of the dynamic occupancy model is that within a primary sampling season the population remains closed, but can be open to local extinction and colonization between seasons (Kéry, 2010). Cusk are described as a "station keeping bottom" or "station keeping cover" species (Auster & Lindholm, 2005), and evidence suggests that the home range of cusk is small (Dultz, 2013).
Furthermore, cusk show strong affinity for complex habitats while avoiding substrata with no structure (Bigelow & Schroeder, 1953;COSEWIC 2003). Cashes Ledge is surrounded on all sides by a number of deep basins, most notably Ammen Basin to the east and Cashes Basin to the west (Figure 2), which consist of unconsolidated substrata (Uchupi, 1966;Uchupi & Bolmer, 2008 Throughout the four sampling seasons within this study, the same camera setup was used. No variability in detection probability should therefore be expected due to gear differences. In investigations using camera surveys combined with fish traps to assess detectability using simultaneous data collection methods (Coggins et al., 2014), the approach of adding cameras to other gear types has been recommended . In the current study, the camera system was used in isolation, and given the limited field of view and a lack of control over orientation once on the seabed, detection probabilities are expected to be and are less than one.
No significant effects were observed for any of the covariates tested in the initial occupancy submodel of model B. While it is likely that some covariates are missing from this submodel, due to the small sample size, it is not possible to add more terms without overfitting the model. The covariates that were included were of importance in two ways. First, they were important in terms of what has previously been reported as driving cusk habitat: depth and surface complexity (Hare et al., 2012). Secondly, maximum current speed and maximum current speed squared were included to test whether hydrodynamics play a role not only in determining cusk behavior, but also in limiting cusk habitat choice. Nevertheless, the main purpose of this study was to give consideration to the detectability issue in marine fish occupancy modeling. Given that the dynamic occupancy model is able to handle constant initial occupancy probabilities across all sites in the study domain, this lack of fit for the initial occupancy state does not present problems for inference about detectability.

| Recommendations for future studies
In this study, detection probability ranged from 0.59 to 0.88 in the noncanopy forming region and from 0.03 to 0.64 in the canopy forming algae region. These results are broadly comparable to detection probabilities from other studies using camera surveys with other modeling approaches to detect marine fish Coggins et al., 2014). These detection probabilities are conditional on cusk being present at the site being observed and also need to be considered in light of the fact that bottom time for the camera was relatively high in this investigation. It should also be noted that, for some species of fish, the presence of camera equipment on the seabed may encourage more curious individuals to investigate so may introduce some bias to the results (Stoner, Ryer, Parker, Auster, & Wakefield, 2008). In both models assessed in this study, the relationship between sampling effort and cusk detection was found to be insignificant. The term was left in the models, however, due to its F I G U R E 6 Tidal curve (top), current speeds (middle) and predicted probability surfaces (bottom) along a full tidal cycle during the summer 2006 sampling season. Surfaces are predicted from the detection model outputs. Numbers of each panel in the gray boxes (t1-t39) correspond to the numbers on the tidal curve. The darker regions in t10 and t30 indicate the areas with canopy forming species of algae based on extinction depth for those species at Cashes Ledge, while the lighter regions are the areas with no canopy forming species of algae theoretical importance; one of the most important factors affecting the detection of any organism is the amount of effort put into trying to detect it. Especially when dealing with small sample sizes, insignificant covariates should be left in models when they are theoretically important (Schuenemeyer & Drew, 2010).
While most of the spatial variability in detection probability can be explained by tidal phase, some small differences persist within each of the bottom types in Figure 6. These can be explained by noninvasive (Bouchet & Meeuwig, 2015;DeCelles, Keiley, Lowery, Calabrese, & Stokesbury, 2017). The need to consider detectability is therefore paramount to obtaining unbiased results.
It is recognized that organisms living in marine environments exist in a multidimensional space, where even though they may live on the seabed, the constant flux of the water column plays an important role in shaping their distributions (Brown, Smith, Lawton, & Anderson, 2011). While the same is true in the terrestrial environment (Jung, Kaiser, Böhm, Nieschulze, & Kalko, 2012), the data needed to describe the nature of the ocean are often harder to collect or generate. Broad scale covariates used in marine investigations, such as temperature, current speed and direction, often come from depth averaged, or surface values of the covariate of concern.
Additionally, if these covariates are not collected at the time of survey as is often the case, they need to be found as records elsewhere, or modeled using an appropriate method. Temperature, for example, is known to affect the metabolism and behavior of marine fish forecasting these types of data may not be a viable option for many researchers planning future studies, it is a recommendation of this study that an effort be made to consider the fine-scale variability of any feature of the environment that may impede detection of their target species.

| CON CLUS ION
This study demonstrates a novel, multifaceted approach to produce a dynamic occupancy model for a species of concern in a closed area. It has generated estimates of occupancy that are considerably greater than occupancy measured from observation alone, for the first time using outputs from a high-resolution 3-dimensional hydrodynamic model in such a modeling framework.
While the need for species distribution models to consider the 3-dimensional nature of the marine environment has been documented previously (Duffy & Chown, 2017), this study reinforces it F I G U R E 7 Observed and median estimated proportions of sites occupied by cusk at Cashes Ledge during the four sampling periods; S1 (Summer 2006), S2 (Spring 2007), S3 (Summer 2007 andS4 (Autumn 2007). Also shown are the 50% and 95% credible intervals for the posterior estimates of the proportion of occupied sites using modern statistical methods. Using the outputs, this investigation has shown how the behavior of cusk changes in two different environments as a function of current speed. This behavior has implications for the detectability of the species, which in turn has implications for the occupancy estimates. This imperative to consider detectability in marine SDM studies is true not only for cusk, but for all species surveyed using noninvasive sampling methods.
It holds especially true for areas where managers must use these methods to monitor stocks. Failing to recognize the limitations of models that do not account for imperfect detection will impact future estimates of abundance, potentially for many species. It is imperative that practitioners of future marine SDM applications consider the detectability of the species under study. In order to do this, they must first understand the processes that govern the fine scale, time-varying features of the environment that may affect detectability, not just for the species in question but also for the specific habitat type being observed  in order to obtain unbiased estimates of occupancy in the marine environment.

ACK N OWLED G M ENTS
This research was part funded by a Fulbright-Marine Institute

DATA ACCE SS I B I LIT Y
Data will be archived on Dryad.