Physics–Guided Advection–Diffusion–Reaction Model for the Environmental Dynamics of Bigeye Tuna in the Eastern Pacific

1. Introduction

Bigeye tuna (Thunnus obesus) are a major pelagic predator in the Pacific Ocean and support one of the world's most valuable fisheries [1,2]. As highly migratory species, their spatiotemporal distribution is exquisitely sensitive to marine environmental gradients and large-scale climate variability [3,4,5]. In particular, the El Niño-Southern Oscillation (ENSO) periodically reshapes the thermal structure, ocean circulation, and primary productivity of the Pacific Ocean, triggering profound reorganizations of pelagic habitats [6,7,8]. For instance, the zonal displacement of the western Pacific warm pool during El Niño events directly drives the large-scale eastward migration of tuna populations [8,9]. Therefore, understanding and accurately predicting the distribution dynamics of bigeye tuna under climate forcing is crucial for sustainable fisheries management and dynamic spatial conservation [10,11,12].

Despite the well-documented impacts of climate on marine ecosystems, accurately simulating tuna distribution remains challenging due to the complex interplay between physiological constraints and the three-dimensional marine environment. Previous habitat studies have primarily focused on sea surface temperature (SST) and dissolved oxygen to define the physiological boundaries of tunas [3,13,14,15]. However, these simplified environmental indicators often fail to capture the habitat compression or expansion experienced by pelagic fish. Cold, hypoxic subsurface layers can compress the vertical habitat available to high-performance tropical pelagic fishes [40]. Bigeye tuna exhibit pronounced deep foraging behaviors, making their distribution highly dependent on the upper ocean's thermal structure, specifically the mixed layer depth (MLD) and ocean heat content (OHC) [4,13,14,16]. The MLD determines the vertical accessibility of the deep scattering layer, while the OHC integrates the thermodynamic state of the water column, acting as a critical energy reservoir during climate anomalies. Relying solely on a limited set of environmental variables overlooks the synergistic constraints of the 3D marine environment, creating a critical limitation in current habitat suitability assessments.

Beyond static environmental suitability, the spatial dynamics of highly migratory species are fundamentally governed by ocean physical transport, yet these processes are frequently isolated from habitat selection in existing models. Population movement is driven by both large-scale advection from ocean currents and effective diffusion representing unresolved mesoscale and submesoscale mixing and dispersal [17,18]. While climate-driven velocity fields force the directed migration of populations, mesoscale eddies create highly heterogeneous diffusion pathways, effectively trapping nutrients and aggregating prey to form localized hotspots [16,19,20]. However, traditional spatial models often treat physical transport and biological habitat preference as independent modules or oversimplify diffusion as uniform random walks. This artificial separation obscures how climate anomalies simultaneously alter passive physical dispersal and active biological aggregation, limiting the ability to mechanistically explain distribution shifts under extreme ENSO events.

Methodologically, current modeling frameworks struggle to reconcile mechanistic physical-ecological processes with the statistical rigor required to assimilate rich observational data. Pure data-driven statistical models (e.g., machine learning) excel at fitting historical abundance patterns but often have limited mechanistic interpretability. They often conflate density-dependent population regulation with environmental adaptation, which can lead to false "niche drift" predictions under future climate scenarios [21,22,23,24]. Conversely, complex mechanistic ecosystem models provide comprehensive physiological descriptions but suffer from limited parameterization capacity and high computational barriers [16,19,20]. There is an urgent need for an intermediate-complexity framework that couples environment-dependent habitat selection with heterogeneous physical transport, while remaining computationally tractable for rigorous data assimilation.

To address these critical gaps, this study proposes a Physics-Guided Advection-Diffusion-Reaction (PhyG-ADR) modeling framework tailored for bigeye tuna in the eastern Pacific. By embedding a multi-factor non-linear habitat selection function into a partial differential equation, the PhyG-ADR model explicitly couples ocean physical transport with density-dependent ecological responses. Furthermore, we leverage Bayesian parameter inference to integrate long-term fisheries observations, directly estimating the key dynamical parameters governing advection, diffusion, and environmental preference. Specifically, this study aims to answer the following core scientific questions: (1) How do upper-ocean environmental structure and physical transport processes synergistically drive the spatial distribution of bigeye tuna? (2) What is the relative contribution of habitat selection versus advection-diffusion during extreme ENSO events? (3) Can the PhyG-ADR framework successfully decouple the species' fundamental habitat preference structure from density-dependent population fluctuations? By answering these questions, this research provides a robust, interpretable tool for assessing the climate adaptability of pelagic fisheries and supporting dynamic ocean management.

2. Materials and Methods

2.1. Data Sources and Preprocessing

We assembled a spatio-temporal database covering the eastern Pacific Ocean (EPO), with a focus on the eastern equatorial region. The study period (1994–2012) spans multiple ENSO warm and cold events and matches the period analyzed by Lian and Gao [4] for comparability. Monthly catch-per-unit effort (CPUE) indices for bigeye tuna were obtained from the Inter-American Tropical Tuna Commission (IATTC)[4]. Only high-quality records with unambiguous spatial and temporal tags were retained. Observations were subjected to standard quality control (removal of duplicates and outliers, unit harmonization, and geolocation checks) and aggregated by month on a 1° × 1° latitude–longitude grid. CPUE reflects relative abundance but is influenced by fishing effort and targeting behavior; aggregation mitigates sampling irregularity, and no spatial interpolation was applied to the fishery observations in order to avoid generating artificial data. The resulting dataset contains approximately 60,000 monthly grid-cell observations of CPUE.

Environmental drivers were assembled from multiple sources to characterize the thermal, dynamic, biogeochemical, and topographic conditions of bigeye tuna habitat. The candidate predictor set included sea-surface temperature (SST), wind speed, mixed-layer depth (MLD), sea-surface height (SSH), chlorophyll-a concentration (CHL-a), ocean heat content (OHC), dissolved oxygen (O2), bathymetry, and marine heatwaves (MHWs) (Fig. 1) . MHWs were identified from daily NOAA OISST v2.1 using a seasonally varying 90th-percentile threshold calculated over the 1982–2011 baseline period. Events lasting at least five consecutive days were retained, and events separated by gaps of no more than two days were merged. Here we used mean intensity to represent the index of HMWs.

SST fields were taken from NOAA’s Optimum Interpolation SST version 2, a global 0.25° daily dataset that merges satellite and in-situ observations. SSH and related surface dynamic conditions were obtained from AVISO/CMEMS multi-mission altimetry. Mixed-layer depth and ocean-heat content, defined as the vertically integrated heat content above the 26 °C isotherm, were taken from NOAA upper-ocean temperature analyses. Dissolved oxygen fields relevant to subsurface habitat conditions were derived from climatological and reanalysis products. Chlorophyll-a was used as a proxy for lower-trophic productivity, bathymetry described large-scale topographic structure, and marine heatwaves were used to characterize extreme thermal anomalies. All environmental variables were averaged to monthly means and regridded to the 1° × 1° grid(Figure 1). To focus on interannual variability, linear trends associated with long-term warming were removed from SST, OHC, and SSH-related fields. We applied identical land-sea masks and carefully treated coastal and boundary cells to avoid artefacts. Unlike the environmental data, CPUE observations were not spatially interpolated; each grid cell retains only real observations.

2.2. Model Formulation

We established a partial differential equation model for the mechanic distribution of tuna species in Pacific Ocean (Figure 2). And we simulated the biomass density C(x,y,t) using an advection–diffusion–reaction equation:

$${\displaystyle \frac{\partial \mathrm{C}}{\partial \mathrm{t}}}=\nabla \cdot (\mathrm{D}\nabla \mathrm{C})-\nabla \cdot \left(\mathrm{uC}\right)+\mathrm{R}(\mathrm{C}, \mathrm{E})$$

(1)

where u(x,y,t) is the advection velocity field, D(x,y,t) is an effective diffusion coefficient, and R(C, E) is a reaction term representing local growth and mortality as a function of biomass and environmental variables E [17,18,25]. The model is solved on the monthly 1° grid using an explicit finite-difference scheme with appropriate boundary conditions (zero-flux at land boundaries and open boundaries at 100°E and 290°E). A monthly time step is used, consistent with the temporal resolution of the observations.

2.2.1. Advection Velocity

Directed movement is represented by an advection term $\nabla \cdot \left(\mathrm{u}\mathrm{C}\right)$. A purely geostrophic velocity would not capture behavioural responses to climate anomalies; we therefore augment the baseline geostrophic field with a climate-driven component estimated from the migration of the population centre of gravity (COG). At each month, the COG of observed biomass is computed as $\mathrm{COG}\left(\mathrm{t}\right)=\left({\stackrel{\mathrm{-}}{\mathrm{x}}}_{\mathrm{t}}-{\stackrel{\mathrm{-}}{\mathrm{y}}}_{\mathrm{t}}\right)$:

$${\stackrel{\mathrm{-}}{\mathrm{x}}}_{\mathrm{t}}={\displaystyle \frac{\sum _{\mathrm{x},\mathrm{y}} {\mathrm{C}}_{\mathrm{aba}}(\mathrm{x},\mathrm{y},\mathrm{t})\cdot \mathrm{x}}{\sum _{\mathrm{x},\mathrm{y}} {\mathrm{C}}_{\mathrm{aba}}(\mathrm{x},\mathrm{y},\mathrm{t})}},{\stackrel{\mathrm{-}}{\mathrm{y}}}_{\mathrm{t}}={\displaystyle \frac{\sum _{\mathrm{x},\mathrm{y}} {\mathrm{C}}_{\mathrm{aba}}(\mathrm{x},\mathrm{y},\mathrm{t})\cdot \mathrm{y}}{\sum _{\mathrm{x},\mathrm{y}} {\mathrm{C}}_{\mathrm{aba}}(\mathrm{x},\mathrm{y},\mathrm{t})}}$$

(2)

The difference between consecutive COGs defines a displacement vector, and the climate-driven velocity ${\mathrm{u}}_{\mathrm{clim}}$is obtained by dividing by the monthly time step $∆\mathrm{t}$. The total velocity is then defined as:

$$\mathrm{u}={\mathrm{u}}_{\mathrm{g}}+{\mathsf{\beta }}_{\mathrm{clim}}{\mathrm{u}}_{\mathrm{clim}}$$

(3)

where ${\mathrm{u}}_{\mathrm{g}}$ is the background geostrophic velocity and ${\mathsf{\beta }}_{\mathrm{clim}}$ is a tunable scaling parameter estimated by Bayesian assimilation. This approach allows the model to follow ENSO-driven zonal shifts evident in the CPUE record.

2.2.2. Diffusion

Small-scale mixing and the aggregated effect of unresolved behaviors are represented by an effective diffusion coefficient:

$$\mathrm{D}(\mathrm{x},\mathrm{y},\mathrm{t})={\mathrm{D}}_0+{\mathrm{D}}_1|\nabla \mathrm{SSH}(\mathrm{x},\mathrm{y},\mathrm{t})|$$

(4)

where ${\mathrm{D}}_0$ is a baseline diffusivity and ${\mathrm{D}}_1$ scales the influence of mesoscale eddies represented by the magnitude of the SSH gradient. The use of eddy-enhanced diffusion reflects observations that tuna aggregations are influenced by fronts and eddies. Both parameters are estimated using the MCMC framework.

2.2.3. Reaction Term and Habitat Suitability

Local population dynamics combine density-dependent growth and environmental carrying capacity. We use a logistic structure in which the intrinsic growth rate r(E) and carrying capacity K(E) depend on environmental drivers and are separated to avoid parameter confounding. The reaction term takes the form

$$\mathrm{R}(\mathrm{C},\mathrm{E})=\mathrm{r}\left(\mathrm{E}\right)\mathrm{C}(1-{\displaystyle \frac{\mathrm{C}}{\mathrm{K}\left(\mathrm{E}\right)}})-\mathsf{\gamma }\mathrm{C}$$

(5)

where γ represents a background loss rate that includes natural mortality and unmodelled processes. Density dependence becomes important as $\mathrm{C}$ approaches K(E), while r(E) controls growth when abundance is low. This formulation is analogous to the Schaefer logistic model, in which the intrinsic growth rate r and carrying capacity K set the production curve and density-dependent term. In the continuous version the density-dependent term appears in the reaction term rather than as a discrete difference[25].

The intrinsic growth rate r(E) reflects the physiological effects of temperature and oxygen. We define:

$$\mathrm{r}\left(\mathrm{E}\right)=\mathsf{\alpha }{[{\mathrm{G}}_{\mathrm{T}}\left(\mathrm{T}\right)·\mathrm{M}({\mathrm{O}}_2\left)\right]}^{{\mathrm{r}}_{\mathrm{p}}}$$

(6)

where α is the maximum growth rate. ${\mathrm{G}}_{\mathrm{T}}\left(\mathrm{T}\right)$ is a Gaussian thermal response centred on the optimum temperature ${\mathrm{T}}_{\mathrm{opt}}$ with width ${\mathsf{\sigma }}_{\mathrm{T}}$, which reflects the unimodal thermal niche of bigeye tuna[13]; metabolic rates and swimming capacity decline when the environment deviates from the optimum. ,$\mathrm{M}\left({\mathrm{O}}_2\right)$ is a Michaelis–Menten saturation function of dissolved oxygen:

$$\mathrm{M}({\mathrm{O}}_2)={\displaystyle \frac{{\mathrm{O}}_2}{{\mathrm{O}}_2+{\mathrm{K}}_{\mathrm{O}}}}$$

(7)

while half-saturation constant ${\mathrm{K}}_{\mathrm{O}}$. Michaelis–Menten functions are widely used to represent biological rates that saturate at high substrate levels. This formulation ensures that growth decreases under hypoxic conditions but approaches α when oxygen is plentiful [13,25].

Habitat suitability S(E) integrates multiple environmental factors to determine the carrying capacity K(E). Following niche theory, we assume that suitability is multiplicative and declines if any factor departs from its preferred range[21,22]. The suitability function comprises four non-linear responses:

$$\mathrm{S}\left(\mathrm{E}\right)={\mathrm{G}}_{\mathrm{T}}\left(\mathrm{T}\right){·\mathrm{L}}_{\mathrm{MLD}}\left(\mathrm{H}\right)·\mathrm{M}({\mathrm{O}}_2){·\mathrm{L}}_{\mathrm{OHC}}\left(\mathrm{Q}\right)$$

(8)

where ${\mathrm{G}}_{\mathrm{T}}\left(\mathrm{T}\right)$ is the same Gaussian thermal response used in r(E), ${·\mathrm{L}}_{\mathrm{MLD}}\left(\mathrm{H}\right)$ and ${·\mathrm{L}}_{\mathrm{OHC}}\left(\mathrm{Q}\right)$ are logistic functions of mixed-layer depth and ocean-heat content, and $\mathrm{M}\left({\mathrm{O}}_2\right)$ is the Michaelis–Menten oxygen term. Logistic functions represent threshold-like responses that saturate beyond certain values. Bigeye tuna tend to prefer shallow mixed layers that facilitate vertical foraging and avoid very deep thermoclines[4,13,14]; logistic functions capture this tendency. Ocean-heat content reflects the integrated thermal energy available in the water column and influences food availability and habitat volume. Using these functions allows the model to allocate carrying capacity to regions with suitable thermal, structural and oxygen conditions. The carrying capacity is then expressed as K(E) = K₀× S(E) × [1 + δ P(E)], where K₀ is a baseline capacity and δ P(E) accounts for productivity enhancements when energy reserves (e.g. primary production) are high.

2.3. Parameter Assimilation via Mcmc

Twelve parameters θ = {D₀, D₁,${\mathsf{\beta }}_{\mathrm{clim}}$, α, δ, γ, ${\mathsf{\sigma }}_{\mathrm{T}}$, ${\mathrm{K}}_0$, ${\mathrm{K}}_{\mathrm{MLD}}$, $K_{\mathrm{OHC}}$, $H_{\mathrm{Opt}}$} control the advection, diffusion and reaction terms. To estimate them and quantify uncertainty we adopt a Bayesian framework[26,27,28]. A likelihood function is defined by comparing simulated biomass C(x,y,t) integrated over each grid cell with the observed abundance ${\mathrm{C}}_{\mathrm{obs}}$ after appropriate scaling:

$$\mathrm{L}(\mathsf{\theta } | \mathrm{Data}) = \prod _{\mathrm{x},\mathrm{y},\mathrm{t}}{\displaystyle \frac{1}{\sqrt{2\mathsf{\pi }{\mathsf{\sigma }}^2}}}\exp (-{\displaystyle \frac{{(\log {\mathrm{C}}_{\mathrm{sim}}(\mathrm{x},\mathrm{y},\mathrm{t},\mathsf{\theta })-\log {\mathrm{C}}_{\mathrm{obs}}(\mathrm{x},\mathrm{y},\mathrm{t}))}^2}{2{\mathsf{\sigma }}^2}})$$

(9)

We assume additive Gaussian observation errors with variance ${\mathsf{\sigma }}^2$. Uniformly informative priors are assigned to each parameter based on physiology and literature values. Posterior samples of $\mathsf{\theta }$ are obtained using the affine-invariant ensemble MCMC (Figure 3) sampler of Goodman and Weare [29], implemented using emcee [30]. This algorithm runs multiple chains in parallel and adapts the proposal distribution to improve sampling efficiency. Each chain consists of 10 000 iterations with a burn-in of 5 000 to allow convergence. Convergence diagnostics include trace plots, effective sample size calculations, and the Gelman–Rubin statistic, which compares variance within and between chains; values near 1 indicate convergence [31]. We also monitor the acceptance ratio and adapt the proposal width to maintain it between 0.2 and 0.5. Parameter posterior means and credible intervals are used in subsequent analyses.

3. Results

The PhyG-ADR model reproduces the observed spatio-temporal variability of bigeye tuna abundance in the eastern Pacific. Hindcasts over 1994–2010 show spatial root-mean-square error reductions of 2.4% compared with benchmark statistical models including the Generalized Additive Model (GAM) [32], Support Vector Regression (SVR) [33], and XGBoost [34] (Table 1). Further performance differences are illustrated using a Taylor Diagram [35] (Figure 4), which demonstrates that PhyG-ADR achieves the highest spatial-temporal correlation and successfully captures the true amplitude of variability (normalized standard deviation near 1.0), enabling skillful prediction of the timing of biomass peaks and troughs associated with ENSO warm and cold phases.

Spatial comparison of the observed and model-predicted mean bigeye tuna CPUE in the eastern Pacific Ocean during 1994–2012(Figure 5). The left column shows the observed mean CPUE fields repeated for row-wise comparison [(a), (d), (g), and (j)]. The middle column presents the corresponding predictions from (b) PhyG-ADR, (e) XGBoost, (h) support vector regression (SVR), and (k) generalized additive model (GAM). The right column shows the prediction residuals, calculated as predicted minus observed CPUE, for (c) PhyG-ADR, (f) XGBoost, (i) SVR, and (l) GAM. Positive residuals indicate model overestimation, whereas negative residuals indicate underestimation. Observed and predicted fields share the same CPUE color scale, and residual maps use a common diverging color scale. RMSE values were calculated over grid cells with valid observations. Land areas are shown in gray.

Consistent with the overall performance metrics reported in Table 1, the spatial comparison further showed that PhyG-ADR provided the closest reconstruction of the observed mean CPUE distribution (Figure 5). All four models reproduced the broad zonal distribution of bigeye tuna CPUE across the eastern Pacific Ocean, including the equatorial high-CPUE band and the large-scale meridional gradient. However, their ability to represent regional-scale spatial variability differed substantially. PhyG-ADR achieved the lowest RMSE (0.0886), followed by XGBoost (0.0908), SVR (0.0993), and GAM (0.1164). Thus, the reduction in RMSE achieved by PhyG-ADR was approximately 2.4% relative to XGBoost, 10.8% relative to SVR, and 23.9% relative to GAM. The relatively small difference between PhyG-ADR and XGBoost indicates that their overall predictive accuracy was similar, although PhyG-ADR produced more spatially localized residuals and better retained the observed zonal structure. The PhyG-ADR residuals were generally centered around zero across most of the study domain, with larger errors confined primarily to localized boundary and coastal regions. XGBoost reproduced the broad spatial pattern but showed more spatially scattered positive and negative residuals. SVR exhibited more extensive negative residuals in the central and southern portions of the study area, indicating systematic underestimation in these regions. GAM generated an overly smooth spatial field and failed to reproduce several finer-scale features apparent in the observations. These results suggest that the process-based constraints incorporated into PhyG-ADR improved the representation of spatial structure, although its numerical advantage over XGBoost remained modest.

Through Bayesian parameter assimilation (Figure 3), posterior estimates indicate that the climate-driven advection coefficient ${\mathsf{\beta }}_{\mathrm{clim}}$ is 0.6 ± 0.1 (95 % credible interval), implying that population movement follows about 60 % of the COG displacement, while the eddy-enhanced diffusion coefficient D₁ is significantly positive, confirming the role of fronts in mixing biomass. The intrinsic growth rate α is estimated at 0.25 month⁻¹ (0.18–0.32 month⁻¹) and the optimum temperature ${\mathrm{T}}_{\mathrm{opt}}$ at 25.5 °C (24.8–26.2 °C), consistent with physiological observations that bigeye tuna thrive between 24 °C and 28 °C. The half-saturation constant for oxygen ${\mathrm{K}}_0$ is 2.8 mL L⁻¹, indicating strong suppression of growth under hypoxic conditions. The posterior distribution of the MLD-response midpoint, H_opt was concentrated around 57 m, indicating that the modeled influence of mixed-layer depth changed most rapidly near this depth. The posterior distribution of K_MLD was comparatively broad, suggesting residual uncertainty in the steepness of the MLD response. The OHC-response coefficient K_OHC was centered near −0.07; under the specified sign convention, this negative coefficient indicates a decreasing modeled response with increasing OHC. These parameters describe associations within the habitat-capacity formulation and should not be interpreted as experimentally determined physiological thresholds (Figure 3).

Mechanistic analysis reveals that habitat selection is the dominant driver of spatial aggregation. When the habitat suitability function S(E) is disabled (by replacing it with a constant), the model underestimates the high-abundance zones and overestimates diffuse distributions. In contrast, removing the climate-driven advection component results in poor representation of east–west shifts during El Niño and La Niña events. During the strong 1997/98 and 2009/10 El Niño events, advection accounted for more than 50% of the model-attributed variance in spatial redistribution, whereas the habitat-related reaction component accounted for more than 70% during ENSO-neutral periods. Sensitivity experiments in which the baseline carrying capacity K₀and productivity scaling δ are varied by ±50 % show that the optimum temperature and oxygen thresholds remain within the posterior credible intervals, demonstrating the robustness of the estimated fundamental thermal niche to changes in carrying capacity (Figure 6). The relative contribution of the modeled processes varied among climate states, with advection becoming more influential during strong El Niño events and habitat-related processes contributing more strongly during ENSO-neutral periods.

At the basin-aggregated scale, the monthly mean relative-abundance index reconstructed by PhyG-ADR broadly followed the temporal evolution of the observed series during 1994–2012 (Figure 7). The model reproduced several major declines and subsequent recoveries, including the pronounced decrease around 2005–2006, indicating that it captured a substantial component of the low-frequency variability in the observations. However, the simulated series generally exhibited a smaller amplitude than the observed series: several short-lived maxima were underestimated, whereas some abrupt minima were only partially reproduced. This temporal smoothing suggests that the monthly 1° process-based model represents broad inter-annual changes more effectively than high-frequency fluctuations, which may additionally reflect short-term fish movement, heterogeneous fishing effort, catchability variation, and observation error. Although some reconstructed changes coincided with major ENSO periods, Figure 7 alone does not establish an ENSO-driven response; such attribution requires an explicit comparison with an ENSO index and event-based statistical evaluation.

4. Discussion

4.1. Physical Interpretability and Ecological Insights

The PhyG-ADR model advances fisheries oceanography by bridging the gap between physical oceanography and ecological modeling in a parsimonious but mechanistically explicit framework. Compared with purely statistical habitat models, it provides robust insights into the relative roles of habitat preference and physical movement in shaping fish distribution, and can therefore be more confidently extrapolated to novel climate states [23,24,36]. The dominance of habitat selection in our posterior analysis is consistent with the view that bigeye tuna maintain a relatively conservative thermal niche, as suggested by growth and metabolic studies [13]. While habitat selection dominates the overall distribution pattern, the inclusion of climate-driven advection is critical to reproduce the large-scale east–west shifts associated with ENSO [4,8,9]. The estimated optimum temperature ($\approx 25.5℃$) and half-saturation oxygen threshold ($\approx 2.8 \mathrm{mL} {\mathrm{L}}^{-1}$) align with laboratory-derived tolerance ranges, supporting the ecological validity of the model [13]. Furthermore, the eddy-enhanced diffusion term underscores the importance of fronts and mesoscale structures in fish aggregation [16,19,20]. This finding highlights the necessity of explicitly representing such physical oceanographic features, as done in PhyG-ADR, to inform dynamic fisheries management tools [10,11,12].

4.2. Parameter Assimilation and Mechanism Robustness

The Bayesian assimilation framework provides a rigorous quantification of parameter uncertainty and addresses identifiability issues that have hampered previous mechanistic models. Convergence diagnostics (such as MCMC trace plots and the Gelman–Rubin statistic, shown in Figure 3) indicate well-behaved posterior sampling [27,28,29,30,31]. However, some parameters, such as the productivity scaling factor $\mathsf{\delta }$, remain weakly constrained in the current setup, as CPUE data primarily reflects relative abundance and provides limited direct information on bottom-up trophic dynamics. Future work could incorporate independent estimates of primary production or prey density to reduce parameter equifinality. Specifically, prey density data would help disentangle the complex effects of temperature and food availability on habitat selection, thereby better constraining $\mathsf{\delta }$ and the baseline carrying capacity ${\mathrm{K}}_0$.

4.3. Implications for Management

As demonstrated by the model's ability to accurately track population shifts during ENSO events, static marine protected areas (MPAs) are insufficient for highly migratory species like bigeye tuna, whose distributions are tightly coupled to dynamic oceanographic conditions [2]. The successful parameterization of climate-driven advection and habitat suitability provides a robust basis for implementing Dynamic Ocean Management (DOM) [10,11,12]. By utilizing near real-time environmental monitoring and the predictive capabilities of the PhyG-ADR model, fisheries managers can dynamically adjust conservation boundaries to protect core habitats while minimizing disruptions to sustainable fishing operations.

4.4. Model Limitations and Implicit Vertical Representation

The use of a 1° grid and monthly time step in this study is a strategic choice to balance computational efficiency with macro-ecological spatial resolution, though it inherently neglects sub-monthly movements and explicit vertical behavior (e.g., diel vertical migration). Nevertheless, the 2D horizontal framework implicitly integrates these processes by incorporating mixed-layer depth (MLD) and ocean heat content (OHC) as proxy covariates in the habitat suitability function. These proxies successfully capture key vertical habitat constraints without requiring full 3D modeling. This limitation is particularly important because ocean deoxygenation and expanding oxygen-minimum zones can reduce the vertical habitat available to pelagic fishes [37,38,39]. The spatially localized residuals of PhyG-ADR, particularly relative to the smoother GAM reconstruction, indicate an improved representation of regional distributional heterogeneity (Figure 5). Extending the model to three dimensions with depth-resolved habitat functions and incorporating fishing effort dynamics would further improve ecological realism. Finally, PhyG-ADR offers a valuable middle ground between purely statistical models and fully coupled, highly complex ecosystem models such as SEAPODYM [16,19,20]. While it significantly reduces the parameter space compared with complex ecosystem counterparts, its mechanistic complexity still increases computational cost. The GPU-accelerated solver and ensemble MCMC used here mitigate this issue, but scalability remains a challenge for global, high-resolution applications.

5. Conclusions

This study developed and evaluated a physics-guided advection–diffusion–reaction model with Bayesian parameter assimilation (PhyG-ADR) to reconstruct the monthly distribution of bigeye tuna in the eastern tropical Pacific during 1994–2012. By integrating environmental suitability, effective transport, diffusion, and density-dependent regulation within a unified partial differential equation framework, the model provides a process-informed alternative to purely empirical approaches. PhyG-ADR achieved the lowest prediction error (RMSE=0.0886), compared with XGBoost (0.0908), SVR (0.0993), and GAM (0.1164). Although its improvement over XGBoost was modest, the model better retained the observed large-scale spatial structure and produced more localized residuals. Bayesian inference yielded ecologically interpretable environmental-response parameters, including an estimated optimum sea surface temperature of approximately 27.26 °C. This estimate should be interpreted as the modeled thermal response of fisheries-dependent relative abundance rather than as a direct measurement of the species’ fundamental physiological niche. The ablation experiments further indicated that removing the habitat-selection component caused the greatest deterioration in the reconstructed spatial aggregation pattern, whereas removing the advective component primarily reduced the model’s ability to reproduce ENSO-associated zonal redistribution. Sensitivity experiments showed that changes in relative carrying capacity substantially affected the magnitude of predicted abundance but had less influence on the normalized shape of the thermal-response curve. These findings suggest that environmental suitability and transport play complementary roles in shaping the mean distribution and climate-related redistribution of bigeye tuna within the model framework.

PhyG-ADR therefore provides a potentially useful framework for dynamic habitat assessment and short-term distribution forecasting. However, the present results are based on fisheries-dependent CPUE and a two-dimensional representation that does not explicitly resolve catchability variation, diel vertical migration, subsurface thermal and oxygen structure, prey availability, fishing behavior, or life-history differences. Accordingly, the inferred processes should not be interpreted as definitive causal mechanisms or direct estimates of absolute biomass. Future work should incorporate standardized fishery indices, tagging or fishery-independent observations, depth-resolved environmental fields, and explicitly modeled fishing effort. Independent spatiotemporal validation, particularly during withheld ENSO and marine heatwave events, will also be necessary before the framework can be applied operationally to dynamic fisheries management and marine spatial planning.