The dataset used to train the network was sourced from simulations produced by the coupled NEMO-BAMHBI modeling system (Grégoire et al., 2008; Grégoire and Soetaert, 2010; Grégoire et al., 2025). NEMO (v4.2) solves the ocean primitive equations and is coupled online to the BiogeochemicAl Model for Hypoxic and Benthic Influenced areas (BAMHBI; Grégoire et al., 2025). This modeling framework simulates the marine carbon, nitrogen, oxygen, and sulfur cycles from bacteria up to mesozooplankton, explicitly representing processes in the anoxic layer through a 28-variable pelagic and 6-variable benthic biogeochemical formulation (Macé et al., 2024; Grégoire et al., 2025).

The coupled model is forced by a regional atmospheric model (MAR) (Gallée and Schayes, 1994; Grailet et al., 2024), which provides high-resolution atmospheric forcing (wind, heat fluxes, precipitation) over the Black Sea. For the hindcast simulation considered here (1950–2014), MAR was forced at its lateral boundaries by global atmospheric fields from MPI-ESM (Giorgetta et al., 2013), ensuring regional meteorological conditions consistent with large-scale climate variability. The resulting atmospheric fields drive the coupled NEMO–BAMHBI system, enabling the simulation of physically and biogeochemically consistent historical conditions in the Black Sea.

The NEMO–BAMHBI configuration used in this study has a horizontal resolution of approximately 0.14^° in latitude and 0.18^° in longitude, corresponding to roughly 15 × 15 km in the Black Sea region. The model outputs are provided at daily temporal resolution and include 59 vertical levels from the surface to 2400 m, with enhanced vertical resolution in the upper ocean to better resolve biologically active layers. This configuration does not include online data assimilation; however, different versions of the NEMO–BAMHBI system have been extensively evaluated against available observations in the Black Sea. In particular, Grégoire et al. (2025) document the validation of BAMHBI simulations against satellite-derived chlorophyll (after 1997), dissolved oxygen, nitrate, and Argo-based chlorophyll observations, supporting the ability of the model to represent key physical and biogeochemical features of the basin.

From the comprehensive set of modeled physical and biogeochemical variables produced by the coupled system, a subset was selected and used as input predictors for the U-Net emulator, as summarized in Table 1. These variables constitute the selected dynamic ocean state fields from the full NEMO–BAMHBI output, while spatial and temporal contextual predictors are introduced separately during data preparation.

The NEMO-BAMHBI dataset was vertically averaged over the upper 50 meters of the water column, a reference depth commonly used as an approximation of the euphotic zone in the Black Sea (Serpetti et al., 2025), where photosynthetic activity and phytoplankton biomass are concentrated. This choice reduces noise from deeper, less biologically relevant layers and ensures that the analysis focuses on surface-driven variability. For each variable, the mean and standard deviation were computed across all spatial grid cells and all time steps in the training periods, and the same normalization parameters were applied to the validation and test sets to avoid information leakage. Other normalization approaches (e.g., log transformations, min–max scaling, and Yeo–Johnson transformations) were explored to better accommodate the diverse distributions of variables, but standardization consistently yielded the most stable and reliable results in preliminary tests and was retained as the default.

Input sequences of length w = 32 days were constructed as model input, with a forecast horizon of q = 16 days. These values were selected based on exploratory experiments evaluating model performance across a range of candidate temporal windows, with details summarized in the Supplementary Material (Section S1.1). Increasing the input window improved forecast skill up to approximately one month, beyond which performance gains saturated and sensitivity to noise increased. Similarly, forecast skill declined rapidly for horizons beyond approximately two to three weeks, consistent with the limited intrinsic predictability of daily chlorophyll variability.

Marginal improvements were occasionally observed for larger input windows when additional mechanisms were introduced to help organize and extract information from the extended temporal context, such as early attention modules (Hu et al., 2018). This suggests that longer input windows may be beneficial when combined with architectures explicitly designed to manage longer-term dependencies, a point that remains open for future investigation. To expand the effective training set while limiting redundancy, input sequences were generated using a sliding window with a stride of 15 days. Smaller strides (e.g., 1 day) led to excessive overlap and overfitting, while fully non-overlapping sequences reduced data availability and degraded performance.

Exploratory feature-selection analyses further guided the final choice of input and output variables”, using a backward feature elimination strategy (Guyon and Elisseeff, 2003). All physical and biogeochemical variables listed in Table 1 were retained as input predictors. While both joint and single-output configurations were explored during model development, the analyses presented here focus on models trained to predict a single target variable at a time (chlorophyll or an individual nutrient). Since chlorophyll is the primary indicator of phytoplankton biomass, the main results and sensitivity analyses focus on chlorophyll for interpretability, with corresponding nutrient results provided in the Supplementary Material.

The dataset was split into distinct periods for training, validation, and testing to ensure robust evaluation. After testing multiple partitioning strategies, the final setup used two non-contiguous training periods (1950–1969 and 1976–2010), with validation data from 1971–1975 and test data from 2011–2014. This configuration ensured that model performance was assessed on data temporally distant from the training set. A one-month buffer period was added between the splits to reduce potential contamination and bias. Although trained on full-basin data, the model was also applied to smaller sub-regions to investigate localized ecosystem variability and drivers. Four pilot sites were selected to capture contrasting physical and biogeochemical regimes within the Black Sea (Figure 2). The northwestern Black Sea shelf is strongly influenced by riverine nutrient inputs, particularly from the Danube River, which enhance coastal productivity and shape local biogeochemical gradients (Ragueneau et al., 2002). Basin-scale circulation in the Black Sea includes western and eastern gyres, a dominant rim current, and mesoscale variability that contribute to regional contrasts in circulation and mixing (Stanev and Beckers, 1999). The Bosporus region reflects exchange dynamics between the Black Sea and the Mediterranean, with distinct effects on stratification and hydrographic structure (Özsoy and Ünlüata, 1997). The southwestern shelf exhibits different shelf–basin interactions compared to the northwestern shelf, while the offshore open-ocean site represents basin-scale conditions where coastal forcing is minimal. These regions reflect contrasting oceanographic regimes representative of major coastal and offshore environments in the Black Sea.

In addition to the physical and biogeochemical input variables, spatial and temporal information was included to help the model learn location-specific and seasonal patterns. Spatial information was provided through input channels for latitude, longitude, bathymetry, and distance to coast. Temporal information was encoded as the day of the year (DOY) using a standard cyclic sine–cosine transformation, day₁ = cos(2π DOY/365) + 1 and day₂ = sin(2π DOY/365) + 1, allowing the network to represent seasonal phase without a discontinuity at the year boundary. This helps the network capture repeating seasonal dynamics across the dataset. The contribution of these spatio-temporal predictors was assessed as part of the backward feature elimination procedure described in the Supplementary Material (Section S1.1).

Model performance is assessed using complementary metrics that quantify both magnitude errors and spatial pattern agreement. Prediction accuracy is evaluated using the root mean squared error (RMSE) and mean absolute error (MAE), which measure the average deviation between predicted and reference values. Spatial pattern correspondence is quantified using the Pearson correlation coefficient, which measures the strength of linear association between predicted and reference maps (Cohen et al., 2009).

To further assess spatial structural similarity beyond pixel-wise (grid-point) errors, the multi-scale structural similarity index (MS-SSIM) (Wang et al., 2003) is employed. MS-SSIM compares predicted and reference fields across multiple spatial scales and accounts for differences in contrast, gradients, and structural organization. Although originally developed for image analysis, MS-SSIM has been used in geoscience and meteorological applications to quantify structural realism and spatial coherence in gridded fields (Harris et al., 2022; Fang et al., 2025).

Forecast skill as a function of lead time is evaluated using a basin-mean spatial RMSE, obtained by averaging the grid-point RMSE over the full basin for each lead day. A relative RMSE is additionally computed by normalizing the RMSE by the mean chlorophyll concentration of the reference field over the evaluated domain (e.g., at each grid cell for spatial maps or over a pilot site for regional time series), yielding a dimensionless, scale-aware error metric. Relative RMSE was preferred over percentage-based metrics such as the mean absolute percentage error (MAPE) or mean absolute percentage deviation (MAPD) because it remains well defined in low-chlorophyll regimes, where normalization by instantaneous values can lead to unstable or misleading errors.

A daily climatology over the period 1950–2014 was constructed as the multi-year mean for each day of the year. To isolate robust seasonal signals, the climatology was smoothed using a Gaussian filter (Steinacker, 2021) with a temporal standard deviation of 3 days and a spatial standard deviation of [2, 2] grid cells. This smoothed climatology serves both as a baseline reference and for the computation of deseasonalized anomalies, which were obtained by subtracting the climatological mean from each daily field prior to model training.

Several persistence-based benchmarks were used to evaluate the emulator relative to progressively more realistic reference forecasts. Pure persistence assumes that the future state remains equal to the most recent observed value at all lead times (Wilks, 2011). An anomaly persistence baseline was also considered, where the initial deseasonalized anomaly is assumed to persist and is added back to the corresponding climatological mean (DelSole et al., 2011; Kirtman et al., 2014). This removes trivial skill linked to the seasonal cycle.

Finally, a damped anomaly persistence model was implemented,

where A(t) denotes the anomaly of the target variable at time t, τ is the forecast lead time, and τ_d is a decorrelation time scale estimated from the temporal autocorrelation of the anomaly time series. This formulation represents the gradual loss of memory in the system and is commonly used as a realistic null model for predictable climate and ocean variability (Newman, 2013; Niraula and Goessling, 2021).

To interpret the relationships learned by the neural network and quantify the influence of environmental drivers on chlorophyll predictions, a combination of global variance-based, derivative-based, and spatial sensitivity analyses was employed. Together, these methods allow not only the identification of the most influential variables, but also the spatial and temporal scales over which they act.

The model demonstrates substantial skill in forecasting chlorophyll changes across the Black Sea (Figure 3). As expected, the prediction error (RMSE) increases with forecast lead time. However, the neural network emulator consistently outperforms all baseline predictions, including climatology and all persistence formulations. After the first two forecast days, the model also surpasses pure persistence, anomaly persistence, and damped anomaly persistence, which assume varying degrees of temporal memory in the initial conditions (Figure 3a). Variability across test samples is illustrated for the emulator to provide context on forecast robustness. Since Climatology is computed by shifting the forecast window in time, it appears constant across lead days and does not reflect increasing forecast uncertainty like the model or persistence do.

During the initial short-range forecast period (< 4 days), the basin-mean RMSE increases approximately linearly with lead time. This is followed by a steeper, weakly nonlinear growth, consistent with the gradual loss of predictability as the model and the real system diverge. At longer lead times, the rate of error growth decreases and the RMSE asymptotically approaches the climatological error level, indicating a transition toward the limit of intrinsic predictability for daily chlorophyll variability. Spatially, the model can capture basin-scale chlorophyll patterns, but with heterogeneous performances (Figure 3b). Lower RMSE values are observed across much of the open ocean and eastern Black Sea, while higher errors are concentrated along the northwestern shelf and western coastal regions. This likely reflects the increased physical and biogeochemical variability in coastal and shelf zones, where shallow depths, strong riverine input (for example from the Danube), and complex circulation generate rapid, localized changes in nutrient supply and phytoplankton blooms. Such dynamic environments are inherently harder to forecast than the more stable, stratified open sea, where variability is lower and thus more predictable.

However, the northwestern coastal regions also exhibit substantially higher absolute chlorophyll concentrations, which further contribute to the elevated absolute RMSE in these areas. When considering the relative RMSE, a different spatial pattern emerges (Figure 3c), with higher relative errors in the transitional zone between the northwestern shelf and the deeper basin, as well as in the dynamically active southern central basin east of the Sinop Peninsula.

While RMSE quantifies absolute error, complementary metrics such as MS-SSIM assess spatial structure preservation, and Pearson correlation captures temporal variability. The model improves RMSE by approximately 39–60% across variables and also achieves gains of +7–12% in MS-SSIM and +11–20% in correlation. This demonstrates that the model not only reduces error magnitudes but also effectively captures coherent spatial patterns and regional dynamics, even when pixel-level discrepancies remain (Table 2). In addition to the basin-scale metrics reported in Table 2, spatiotemporal validation results for nitrate and ammonium are presented in the Supplementary Material (Section S1.3).

Among the nutrient variables, phosphate exhibits the strongest relative improvement across all metrics, with a 58.8% reduction in RMSE and a 20.5% increase in Pearson correlation compared to climatology (Table 2). In contrast to chlorophyll, whose RMSE increases rapidly beyond the first few forecast days, phosphate errors increase approximately linearly and more gradually with forecast horizon. Although mean forecast skill for nitrate and ammonium remains high, both variables exhibit a wider spread in prediction variability across test samples, indicating reduced forecast stability relative to phosphate (Supplementary Figures S1 and S2). These differences may suggest that phosphate follows more stable and slowly evolving dynamics, while nitrate and ammonium are more sensitive to episodic and spatially heterogeneous processes.

It is worth noticing that large areas of the basin maintain persistently low chlorophyll levels, naturally limiting RMSE and SSIM sensitivity when averaged over the full basin performance while significant improvements relative to the considered baselines are particularly pronounced in biologically active regions. These results confirm that the model delivers meaningful improvements over baseline methods, producing forecasts that are both more accurate and more ecologically realistic at sub-basin and basin scales.

Overall, the results highlight that the emulator is able to minimize absolute errors relative to the outputs of the process-based model, while maintaining both structural and dynamic fidelity in its forecasts. These are critical attributes that support the applications of this emulator in short-term predictions and driver analyses, something which can contribute to more informed and effective ecosystem management in the Black Sea.

All basin-scale and spatially averaged results reported here are based on the 2011–2014 test period and on sensitivity analyses computed for the emulator trained with the final selected configuration and temporal split, providing a consistent reference framework for interpreting driver importance. In addition to forecasting chlorophyll and nutrient distributions, the emulator model enables analyses of the key environmental drivers influencing these predictions. Specifically, it enables the assessment of which features most directly influence emulator predictions, either alone or in combination with other inputs. Sobol indices have been calculated across the entire Black Sea, as well as within each of the four pilot sites to investigate regional and seasonal differences in dominant drivers. In this context, sensitivity metrics are interpreted at the level of individual predictors, acknowledging that some inputs are physically co-varying, in order to facilitate interpretation of emulator behavior (see Supplementary Section S1.2). Alternative approaches, such as grouping correlated variables prior to sensitivity analysis, could provide complementary perspectives, but would reduce variable-level interpretability, which is a central objective of the present analysis.

The sensitivity analysis demonstrates that chlorophyll is the dominant driver for future chlorophyll predictions that typically explains > 50% of the variance of the model across the basin and throughout the year (Figure 6). This is not surprising and indicates strong temporal autocorrelation in the dataset. The influence of chlorophyll is highest for short time predictions (< 2 days) and decreases toward q = 16 days. This decline in feature importance is even more pronounced for the secondary drivers like PAR, whose contribution nearly halves from the first to the last timestep of the forecast horizon. Conversely, important drivers such as ammonium and phosphate show the opposite pattern. Indeed they have negligible impact at the initial timestep but explain 10% of total variance for the last forecast day. Physical variables such as temperature, salinity, mixed layer depth, and shortwave radiation have only modest contributions (1–2%), while other drivers, including ocean currents, SSH, and vertical velocity, are negligible (< 1%).

A neural network-based emulator (i.e., surrogate model) has been developed and tested to forecast the output of a process-based ocean–biogeochemical model for chlorophyll and nutrient concentrations in the Black Sea. The results show that the model can reproduce the spatiotemporal variability of both chlorophyll and nutrients with high accuracy and low errors. This suggests that the machine learning model has effectively learned the key ecosystem dynamics simulated by the more complex process-based ocean model for the Black Sea.

The availability of an accurate emulator offers several key advantages in predicting and understanding dynamics in the Black Sea. First, it dramatically reduces computational costs, enabling faster simulations that are particularly important for real-time forecasting and scenario analysis. For reference, the NEMO–BAMHBI configuration used in this study (15 km horizontal resolution, without data assimilation) requires on the order of one hour of wall-clock time to simulate one model year using approximately 192 CPU cores. In contrast, once trained, the neural network emulator produces a 16-day chlorophyll forecast in a fraction of a second on a single GPU, while training the emulator requires on the order of a few hours on a modern GPU. This represents several orders of magnitude reduction in computational cost for inference compared to the original process-based model. Second, it facilitates extensive sensitivity analyses and ensemble simulations that would be very expensive with the original process-based model. Sensitivity analyses including variance-based (Sobol), derivative-based (DGSM), and elasticity-based locality metrics can be used in combination to analyze the feature importance of the input variables used in the model, thus helping to understand ecological surface-layer drivers of specific ecosystem components. Additionally, the emulator can capture essential dynamics of the marine ecosystem, hence it can serve as an effective tool for exploring responses of the system to changing environmental conditions or management interventions. This is particularly important for supporting decision-making in marine resource management and climate adaptation planning.

The developed emulator, trained to forecast multivariate data up to q = 16 days ahead, provides considerable skill for short-term forecasting and adaptive environmental management in the Black Sea. The 16-day forecast window considered here is illustrative in that respect. Indeed, the emulator framework can be tailored to specific regional or application needs, which can be done by increasing the length of the input window w, adjusting the forecast horizon q, or moving from chlorophyll and nutrients to other target variables. On the other hand, the suitability of a given forecasting range also depends on the local temporal dynamics and predictability. In more stable regions, longer forecasts may remain reliable, while highly dynamic areas may benefit from shorter high-frequency updates. The use in inputs of high-resolution observations can also be considered to improve the characterization of local dynamics and to improve predictive skill. This theme has not been directly addressed in this paper, but it may provide an interesting future research direction.

Within the flexible framework, the emulator can support a range of potential, model-based application domains. When trained on oxygen fields from a biogeochemical model, short-range forecasts could be explored as part of numerical experiments or prototype early-warning systems for hypoxic conditions in the Black Sea. Similarly, emulator-based forecasts of modeled nutrient dynamics (e.g., nitrate and phosphate) may be used to investigate the sensitivity of nutrient distributions to episodic inputs, such as rainfall-driven runoff or seasonal agricultural activity, within scenario-based analyses. In fisheries and ecosystem management contexts, the emulator may serve as a rapid exploratory tool to assess how modeled environmental conditions relevant to spawning habitats or migration corridors respond to changing physical or biogeochemical drivers.

A promising application of this framework may lie in the investigation of extreme or rapidly evolving events, such as harmful algal blooms (HABs), which pose serious threats to marine ecosystems, fisheries, and coastal livelihoods. However, the present study does not explicitly evaluate predictive skill during strongly anomalous years, and the temporal test period primarily reflects typical variability rather than rare extremes. In addition, part of the short-term forecast skill arises from the strong autocorrelation of chlorophyll and nutrient fields, which may reduce predictive capability during abrupt regime shifts or externally forced events.

Within these limitations, the emulator could nevertheless be used to explore the emergence of bloom-like conditions in numerical simulations by identifying modeled precursors, such as rapid chlorophyll accumulation or elevated nutrient availability. Recent work has demonstrated the feasibility of ML-based HAB prediction. For instance, Gupta et al. (2023) developed a machine learning model to forecast cyanobacterial bloom intensity in Lake Erie at sub-monthly lead times using meteorological and remote sensing data. On a broader scale, Wu et al. (2025) proposed a global deep learning framework to predict surface chlorophyll concentrations from physical oceanographic inputs, demonstrating strong performance across multiple marine regions. Within this context, the emulator developed here provides a complementary basin-scale forecasting capability for both chlorophyll and nutrient dynamics, enabling future HAB early-warning applications.

The emulator was trained on NEMO–BAMHBI simulation data spanning 1950–2014, which naturally constrains what it can learn about the behavior of the system. This ∼65-year period captures a range of historical conditions, including the mid-20th century eutrophication and subsequent ecosystem regime shifts in the Black Sea (Mee, 1992; Oguz and Velikova, 2010).

It must be noted that the emulator is not simply interpolating patterns from the past but learns underlying statistical and dynamical relationships between key biogeochemical and physical variables as represented in the NEMO–BAMHBI model. This enables the emulator to capture consistent dependencies and interactions within the NEMO–BAMHBI system, allowing it to generate plausible responses under a range of input conditions that remain within the domain of the training data. As a result, the emulator can be used to explore and anticipate potential future states of the Black Sea ecosystem within the modeled framework, particularly when these remain similar to the historical variability represented in the training period. However, this predictive capability is inherently constrained by the boundaries of the training data (1950–2014). While the emulator may generalize to near-future scenarios, its reliability diminishes when exposed to forcing conditions or extreme events that lie well outside the historical range. Hence, while the emulator provides a powerful tool for efficient scenario testing and short- to medium-term forecasting, caution must be exercised when interpreting results under strongly non-stationary future climates.

To address these limitations, future work could incorporate simulation outputs extending to future scenarios, comparing model performance when trained on different time periods. Training the emulator on both historical and projected data would allow it to learn evolving ecosystem responses under climate change as represented in numerical models and improve its robustness when applied to novel conditions. This would also enable investigation of how dominant phytoplankton drivers, and their relative importance as inferred from sensitivity analyses, may shift over time, for instance from nutrient-driven dynamics in the past to increased control by temperature, light, or salinity in a more stratified and saline future (Pörtner et al., 2014; Podymov et al., 2021). An additional important extension would be to explicitly enforce biogeochemical mass-balance constraints within the neural network, for example by ensuring consistency between phytoplankton growth and nutrient uptake according to Redfield-type stoichiometric relationships (Redfield, 1958; Geider and La Roche, 2002). Such physics- and ecology-informed constraints could further enhance the physical realism and long-term stability of emulator-based forecasts and could be implemented within a physics-informed neural network (PINN) framework (Raissi et al., 2019; Karniadakis et al., 2021; Battina et al., 2025).

The spatial and vertical resolution of the surrogate model is inherently constrained by the resolution of the underlying NEMO–BAMHBI simulation data, which provides daily outputs at approximately 15 km horizontal resolution. While sufficient for capturing mesoscale and large submesoscale variability, this resolution limits the ability of the model to represent fine-scale coastal processes such as narrow boundary currents, sharp frontal gradients, and localized river plumes. Such sub-grid features may be smoothed or underrepresented, potentially leading to an underestimation of extreme values and ecological hotspots. Vertically, this study averaged the upper 50 m of the water column to focus on the photosynthetically active layer and reduce model complexity. While this approach emphasizes the biologically productive zone, it inevitably obscures vertical structure. The Black Sea is known for exhibiting a pronounced deep chlorophyll maximum (DCM) around 30–40 m depth, which can contribute substantially to the total chlorophyll inventory (Ricour et al., 2021). Averaging across this depth range may flatten key features associated with stratification and vertical nutrient gradients, reducing model sensitivity to processes such as nutrient entrainment from below the mixed layer. Alternative strategies such as mixed-layer-informed averaging or retaining vertical resolution in a simplified multi-layer format could offer improved representation of subsurface dynamics without compromising computational efficiency. Future work may explore these directions to further enhance emulator accuracy and ecological realism at the basin scale. Consequently, interpretations of sensitivity and feature-importance results should be understood as reflecting correlations within the vertically averaged surface layer and their consistency with known physical and biogeochemical processes, rather than as a direct resolution of underlying vertical mechanisms such as mixing or nutrient entrainment.

Beyond resolution-related constraints, an additional limitation concerns the choice of input predictors. In addition to physical and biogeochemical predictors, the emulator was provided with explicit spatial (latitude, longitude, bathymetry, distance to coast) and temporal (cyclic day-of-year) context. These variables were evaluated implicitly during backward feature elimination and were retained because their removal resulted in a modest but systematic degradation of validation performance, particularly affecting spatial coherence and seasonal phase alignment. While a 3D U-Net can infer relative spatio-temporal structure from multivariate input fields alone, explicit context variables reduce ambiguity in a geographically heterogeneous basin such as the Black Sea. Importantly, these predictors act as conditioning information rather than mechanistic drivers, and their inclusion may reduce generalization under strongly nonstationary or climate-change conditions where spatial–seasonal relationships evolve beyond the training distribution. For this reason, sensitivity analyses in this study emphasize physical and biogeochemical variables, and future work will explicitly isolate, quantify, and compare the contribution of spatio-temporal context predictors, both individually and in grouped form, to better assess their influence, interpretability, and potential spatio-temporal biases.

A related consideration is that several candidate predictors are partially redundant or physically co-varying, which can complicate attribution and mask more physically interpretable drivers. As mentioned in Supplementary Material (Section S1.1), most available outputs from NEMO–BAMHBI were initially included as potential inputs. A stepwise elimination approach was then used to remove variables that did not improve model performance, retaining only those that showed a measurable effect. While this strategy helped identify the most influential drivers, it also introduced the risk that redundant inputs could mask the importance of more physically meaningful ones.

For example, SWR was initially retained due to its moderate importance scores. However, SWR is not a direct forcing term in the chlorophyll dynamics; its main role is as a surface heat source in the temperature equation of NEMO–BAMHBI. Since temperature was already included as an input, SWR may have partially masked the relative importance of temperature in the model. Furthermore, PAR captures the biologically relevant component of incoming light more directly than SWR. This suggests that the sensitivity attributed to SWR primarily reflects shared variance with thermodynamic and radiative predictors rather than an independent biological effect.

Moreover, the sensitivity analyses, including variance-based Sobol indices, derivative-based DGSM, and elasticity-based locality metrics, revealed that the importance of different drivers varies markedly across regions and seasons. For instance, nutrient variables such as ammonium, nitrate, and phosphate exerted greater influence at later forecast lead times, particularly in coastal areas and during spring bloom periods, whereas light-related variables such as PAR dominated during summer. This spatiotemporal heterogeneity underscores the value of regionalized model interpretation: chlorophyll dynamics in the Black Sea are not governed by uniform processes, but by distinct environmental drivers that vary with both location and time. Recognizing this variation is important for designing targeted monitoring strategies and region-specific management interventions.

Recent advances in hybrid machine learning approaches have shown promising results for improving the transferability of surrogate models to real-world observations. Neural networks pre-trained on large numerical model simulations and subsequently fine-tuned on observational data, such as satellite-derived sea surface height, SST, chlorophyll, or in situ measurements, have been shown to achieve superior predictive skill, improved generalization, and greater physical consistency than models trained solely on either data source (Buongiorno Nardelli et al., 2022; Shin et al., 2024; Martin et al., 2024; Wu et al., 2025). For instance, hybrid training has been successfully applied to sharpen altimeter sea level maps (Archambault et al., 2024), upscale air–sea CO₂ fluxes (Kim et al., 2024), and reconstruct subsurface properties (Lentz et al., 2025). In the study by Amadio et al. (2024), a neural network was used to reconstruct nitrate profiles from Argo oxygen data, and the output was assimilated into a Mediterranean biogeochemical model, improving basin-scale nitrate and oxygen dynamics.

Applying a similar approach in this study, by adapting the emulator architecture to align better with available observational data and fine-tuning it on real chlorophyll or nutrient measurements (such as satellite ocean color or Argo profiles), may reduce simulation-derived biases and enable forecasting directly from observation-based inputs, thereby bridging the gap between model training and operational use. However, such posterior neural network–based assimilation should not be seen as a replacement for established data assimilation systems like 4DVAR (Carrassi et al., 2018), which have a better theoretical basis and are physically rigorous. Its main advantage lies in practicality, as it allows additional observations to be incorporated without requiring access to the source code of the numerical model, adjoint system, or involvement of model developers, which is often time-consuming and costly.

Finally, as computational power and GPU capabilities continue to advance, modeling the entire water column, or at least a part of it, using neural networks may become feasible. This approach would allow for a more comprehensive representation of vertical ecosystem dynamics, eliminating the need for depth averaging and potentially improving the ability of the surrogate model to simulate processes occurring throughout the full depth range of the Black Sea. Exploring these possibilities in future studies could further enhance the robustness and applicability of neural network-based approaches for marine ecosystem modeling.