Supervised and unsupervised machine learning methods for modeling current and future habitat of Peruvian anchovy

Understanding the drivers and potential impacts of environmental variability on the distribution of Peruvian anchovies, the largest single-species fishery on the planet, is essential for their proper management in a changing world. However, the intricate interactions of these organisms and environmental variability require the use of complex models such as machine learning methods. In this study, we compared three methods for producing habitat maps of anchovies: the traditional Generalized Additive Models, the XGBoost, which is a form of supervised machine learning, and a new method based on clustering water types as a form of unsupervised machine learning. We optimized the three methods with a parameter grid search algorithm and compared their capability to replicate the mean state of anchovies by comparing them with presence-absence observations along the Peruvian coastline between 1990 and 2010. We used the output of a physical-biogeochemical model as input for the habitat models to produce distribution maps of anchovy. All models successfully simulated the distribution of anchovies along the Peruvian coastline. Anchovies distribution was concentrated in the southern part of the region during year 1998, characterized by the impact of a canonical El Niño, while their distribution extended until the north of Peru and south of Ecuador in the other years of the study period. We then applied the models to predict potential changes in the distribution of anchovies under projected temperature and wind patterns by the end of the century. We observed a reduction in the probability of anchovy occurrence under conditions of higher temperature and weaker winds. Two of the three habitat models predicted a severe maximum decline by 75% (GAM) and 60% (XGBoost) whereby the clustering model predicted a smaller maximum decline in anchovy occurrence by less than 15%.

Understanding the drivers and potential impacts of environmental variability on the Peruvian anchovy (Engraulis ringens) is crucial to assess the potential climate change effects on this fishery of global importance as well as the economic consequences. The Peruvian anchovy is the largest single-species fishery in the world (Chavez et al., 2003, 2008), accounting for 24 (23 for fish oil) to 47% of the global fishmeal and fish oil supply (Fréon et al., 2017; Avad́ı et al., 2014) which are an important feed source for aquacultural production (Fréon et al., 2014). This species inhabits the coastal waters within the Northern Humboldt Current System (NHCS), typically spanning latitudes from about 4°S up to 24°S in two well-defined stocks (Bouchon et al., 2021). Anchovies respond in complex ways to the strong environmental interannual variability of the southeastern Pacific Ocean. The anchovy fishery was a proliferating business that increased over the 1960's reaching landings of 10.9 Mt in 1970 and collapsing shortly afterwards during the strong El Niño event 72/73 (Alheit and Niquen, 2004). Fish populations declined by over 50% during El Niño 82/83 event (Barber and Chavez, 1983) and catches dropped to only 1 Mt in 1998 (Ñiquen and Bouchon, 2004), possibly due to a re-distribution of the fish toward the more coastal waters due to El Niño 97/98 (Bertrand et al., 2004).

Long-term trends in environmental factors due to climate change might substantially impact the habitat of anchovies. Under climate change, the upwelling center is projected to shift poleward (Rykaczewski et al., 2015), potentially resulting in a displacement of the anchovy population (Pinsky et al., 2013; Checkley et al., 2017). However, any wind changes in the NHCS under climate change and their effects on upwelling strength and coastal sea-surface temperature are still highly uncertain, as the effects of increased stratification due to surface warming and increased upwelling-favorable winds may balance each other to some extent (Gutiérrez et al., 2011; Echevin et al., 2012, 2020). Furthermore, climate change is expected to increase the frequency and intensity of El Niño events (Gergis and Fowler, 2009; Shin et al., 2022) which have affected the fishery in the past.

The drivers of anchovy distribution and its response to environmental changes are not yet completely understood. Chavez et al. (2003) concluded that anchovy abundance may be correlated with cold water regimes, while Bertrand et al. (2011) point out that oxygen also plays an important role. Brochier et al. (2013) and Flores-Valiente et al. (2023) suggest that early stages of anchovy may benefit from climate change as a result of increased retention and potentially higher food availability due to a strengthened stratification, nonetheless, Checkley et al. (2017) note that a future reduction in nutrient supply could diminish plankton productivity and alter its composition, thus, as a source of food for anchovies, impacting their population. Currently, anchovies have a competitive advantage compared to other species due to their ability to cope with low oxygen concentrations and have, therefore, been thriving in decades characterized by such conditions (Bertrand et al., 2011). However, under future environmental scenarios with even lower oxygen anchovies might be replaced by smaller fish (Salvatteci et al., 2022). This indicates a complex interaction between anchovy and environmental dynamics in the NHCS.

Habitat models, also known as species distribution models or habitat niche models, predict potential areas of distribution of marine organisms, such as anchovies. By finding relationships between their observed occurrence and environmental drivers, these models offer valuable insight into potential habitat shifts driven by climate change. To do so, habitat models can be trained using several environmental predictors such as temperature, bathymetry and wave energy, among others (e.g., Luján Paredes, 2016; Oliveros-Ramos, 2014; Gogina et al., 2016; Schubert et al., 2015; Sequeira et al., 2014) to predict biomass, abundance, presence only and presence-absence data (Guisan and Zimmermann, 2000; Oliveros-Ramos and Shin, 2023; Sequeira et al., 2014; Grüss et al., 2014). The relationship between these various predictors and the target variables can be very complex and non-linear, requiring machine learning or other advanced modeling techniques to be captured. Generalized Linear Models (GLM) and Generalized Additive Models (GAM) are commonly used methods for modeling habitat due to their capability of fitting non-normally distributed data and capturing the joint impact of several predictors (Guisan and Zimmermann, 2000). They have been applied in several marine ecosystems (Quinci et al., 2022; Gordó-Vilaseca et al., 2021; Muhling et al., 2020; Murase et al., 2009) including the NHCS (Oliveros-Ramos, 2014).

In recent years supervised machine learning algorithms such as random forests and boosted regression trees have become popular methods with the rise of artificial intelligence applications to scientific problems, as an alternative to the traditional statistical models. The Maxent model uses a maximum entropy approach to model species distribution from presence-only datasets (Phillips et al., 2006). It is one of the most commonly used tools in ecology (e.g., Silva et al., 2018) and, in contrast to R and Python packages, it is distributed as a stand-alone software with a graphical interface (https://biodiversityinformatics.amnh.org/open_source/maxent). Boosted regression trees are another common method in habitat niche modeling (Quinci et al., 2022; Muhling et al., 2020) due to their high predictive performance and interpretability (Elith et al., 2008). The XGBoost is a scalable machine learning system for tree boosting (Chen and Guestrin, 2016), which implements randomization techniques to reduce overfitting (Bentéjac et al., 2021). It has been shown to have a better performance for modeling species distribution compared to other machine learning algorithms (Effrosynidis et al., 2020). Clustering is a method of unsupervised machine learning in which patterns in the data are identified. Unsupervised learning is commonly used for exploratory data analysis but, to our knowledge, has not been applied for making fish habitat predictions so far. Machine learning methods are suitable to disentangle the complex and non-linear links between environmental parameters and to predict the distribution of marine organisms; hence they might be able to capture the intricate relationship between anchovy habitat preferences and its drivers.

The aim of this study is to understand how environmental variability affects the habitat of Peruvian anchovy in the Northern Humboldt Current System and how the probability of finding anchovies might change under climate change scenarios using different habitat models. To do so, we compared three methods; we applied XGBoost and GAM, both supervised learning methods, and a new approach with unsupervised learning that involves clustering. In-situ observations were used to train the models. For climate change estimates, we used a physical-biogeochemical model to predict the anchovy habitat. The model parameter validations were done with model data, hence accounting for biases that might arise from the differences between the physical-biogeochemical model output and field observations. We designed a common algorithm for training and comparing the three methods. We compared the predictions by our three types of models during normal years between 1990 and 2010 and 1998 when anchovy catches were low due to El Niño (Bertrand et al., 2004). Finally, we used our models to simulate potential anchovy distribution under idealized climate change conditions.

Data description

The observational data was obtained from acoustic surveys performed by the Instituto del Mar del Perú (IMARPE) along the Peruvian coastline in transects up to 200 nautical miles (nm) long from the coast toward the open ocean (Figure 1) from 1990 to 2010 (Gutiérrez et al., 2012). The cleaned dataset consists of 219 thousand samples. Anchovy abundance data is highly skewed to the absence class with values as high as 182.5 thousand but a mean of only 216 m 2 /square nautical mile (nm 2). Therefore, we transformed anchovy biomass to presence-absence by assigning the value of 1 (presence) to all non-zero data points. This resulted in 25.3% of the samples with anchovy presence and 74.7% with absence (Figure 1). We transformed all environmental variables that we would later use to train the habitat models with supervised and unsupervised methods to a range from 0 to 1 using the StandardScaler function of scikit-learn (Pedregosa et al., 2011). In addition, for the supervised learning models, we transformed the distributions of the environmental variables with a quantile transformer into a Gaussian-like distribution (Pedregosa et al., 2011).

Ocean physical-biogeochemical model

We employed a coupled ocean physical-biogeochemical model to investigate how the environmental conditions affect the distribution of anchovy. The model combines the CROCO (Coastal and Regional Ocean COmmunity model, Shchepetkin and McWilliams, 2005) for physical environmental conditions with the BioEBUS (Biogeochemical model for the Eastern Boundary Upwelling Systems, Gutknecht et al., 2013) model. The CROCO framework provides the background for simulating ocean circulation, temperature, and salinity on a regional scale that extends from 10°N to 33°S and from 69°W to 118°W with a horizontal resolution of 1 12 ° and 32 vertical sigma levels. However, for this study we used the model output from a subregion covering our area of interest along the Peruvian coastline from 3°S to 19°S and from 70°W to 84°W (Figure 1). The BioEBUS model is integrated with CROCO to represent the biogeochemical environment with a complex nitrogen cycle including oxygen-dependant processes like nitrification and remineralisation. The physical component of the model runs in time-steps of 500 seconds, the biogeochemical component in steps of 1,500 seconds and the model output has a monthly resolution. Here, to align with the acoustic cruise data detailed in Section 2.1, we have narrowed the focus of our study to a smaller domain covered by the available observations, ranging from 3°S to 19°S latitude and from 70°W to 84°W longitude, as illustrated in Figure 1.

Hindcast simulation

To simulate the historical interannual variability in the environment that potentially affects fish distribution, we used an interannual configuration of CROCO-BioEBUS. Surface heat and freshwater flux forcing for CROCO are provided by CFSR (Climate Forecast System Reanalysis) data, and the wind forcing by the CCMP (Cross-Calibrated Multi-Platform) product with a temporal resolution of 6 h. Initial and boundary conditions are from monthly SODA (Simple Ocean Data Assimilation, Carton et al., 2018) and monthly climatology CARS (CSIRO - Commonwealth Scientific and Industrial Research Organization Atlas of Regional Seas, Ridgway et al., 2002). A 30-year spin-up period, by repeating the forcing of the year 1990, ensures model equilibrium before conducting simulations for the period from 1990 to 2010. After spin-up, the model is forced by monthly interannually-varying forcing during the period from 1990 to 2010, which was later used for training the habitat models. The same model set-up was previously used in an End-to-End model study investigating the bottom-up impact of plankton dynamics on higher trophic levels, including anchovy (Hill Cruz et al., 2022). Model evaluations are available in (Hill Cruz et al., 2022; Xue et al., 2024a). Model output is available in (Xue et al., 2023).

We sampled the model output at the same spatial and temporal coordinates where field observations were collected. Sea surface temperature (SST), distance from the coast and depth are well correlated between the model and observations (Table 1; Figure 2). Meanwhile, sea surface salinity and oxycline depth have a lower correlation of 0.63 and 0.56, respectively. Furthermore, chlorophyll a has higher values in CROCO-BioEBUS than in the observations (Appendix 1) and it exhibits a correlation of only 0.28. Comparing the natural logarithm of CROCO-BioEBUS and observations chlorophyll a increased the correlation to 0.4 (Table 1); however, since this is still very low, we did not include it in our habitat models (see Section 2.3). We also did not include distance from the shelf in the habitat models due to its high correlation with distance from the coast (0.88) and depth (0.81; Figure 3). A bi-modal distribution is observed in SST (Figure 2) due to the sampling taking place in two main seasons (see Appendix 1). On the other hand, the two peaks are less distinct in the CROCO-BioEBUS model SST, which also shows a small systematic warm bias relative to observations, consistent with Xue et al. (2024a). One process that could contribute to this warm bias is the absence of chlorophyll-dependent light attenuation in the model: without shading, solar radiation penetrates deeper, warming source waters later drawn to the surface by coastal upwelling (Echevin et al., 2021).

Climate change experiments

We performed a set of three sensitivity experiments (TEMP+, WIND+, and WIND-) to test how our predicted anchovy distribution responds to different environmental states simulated by the coupled physical-biogeochemical model. The sensitivity experiments were designed to represent wind speed changes and warming which could realistically occur by the end of the century according to climate models. These changes were applied to a climatological reference simulation (REFERENCE), which was set up by averaging the forcing of the hindcast simulation described above between 1990 and 2010. The reference simulation was run for a 20-year spin-up period, after which the sensitivity experiments were started. The reference simulation and the sensitivity experiments were then run for a period of 10 years and the last year was used for analysis. Since the climatological forcing contains no interannual variability, the model already reaches a quasi-steady-state within this time-frame in the upper 200–300 m of the ocean, which are relevant for the anchovy habitat.

The warming experiment (TEMP+) was designed as an idealized sensitivity study to understand the mechanistic impact of temperature change on anchovy habitat, reflecting an increase in surface air temperature and water column stratification consistent with end-of-century projections from climate models, while keeping the wind forcing constant. We applied a spatially uniform ocean temperature offset at the open boundaries and to the initial temperature fields, in the form of a vertical profile. To compute this average profile of ocean warming, we took the temperature difference between 1990–2010 as the baseline period and 2080–2100 as the end-of-century target from the CMIP5 (Coupled Model Intercomparison Project Phase 5, Taylor et al., 2012) multi-model ensemble average under the RCP8.5 emission scenario and averaged it spatially over our ocean model domain. The 1990–2010 baseline was chosen to match the period covered by both the observational dataset and the CROCO-BioEBUS hindcast simulation. We then applied the surface value of the previously computed water temperature change profile (2.2 °C) as a spatially uniform offset to the surface air temperature, which is used by the model for computing air-sea heat and momentum fluxes. As this experiment aims to probe the mechanistic response of habitat suitability to warming rather than to deliver a climate change projection, the precise choice of emission scenario is secondary; what matters is that the perturbation magnitude falls within the range of projected warming for this region. A similar approach of modifying surface air temperature using the bulk forcing option in CROCO has been used before in regional climate change downscaling simulations to add a forcing perturbation (e.g., Espinoza-Morriberon et al., 2017; Echevin et al., 2020).

CMIP models typically have a coarser resolution (about 1°) than the regional CROCO BioEBUS configuration used here (1 12 °), and, therefore, tend to fail to capture the fine scale characteristics of the coastal upwelling region (Small et al., 2015). By increasing air temperature according to CMIP projections, itself designed as an idealized sensitivity study, we allowed the regional model to adapt to the change in forcing and develop its own circulation response. Consequently, despite the spatially uniform increase in air temperature, the resulting sea surface temperature change is not uniform across the model domain but also depends on the physical response of the NHCS. Additionally, because the ocean model does not feedback to the atmosphere, directly imposing a change in sea surface temperature, rather than air temperature, would create a heat imbalance between ocean and atmosphere, producing a persistent, unrealistic heat flux out of the ocean.

To further investigate the uncertainty of potential wind changes and their effect on the anchovy habitat in the region, we increased the wind speed by 10% (WIND+) in a first sensitivity simulation and decreased it by 10% (WIND−) in a second one. This corresponds to a wind stress change of approximately ± 20%, as the wind speed enters quadratically into the wind stress parameterization. All surface fluxes (momentum, heat, CO 2, and O 2) are calculated by a bulk formula and are thus consistent with these wind changes. Note also that the wind speed changes are small compared to some previously published sensitivity studies (e.g., Echevin et al., 2012; Travers-Trolet et al., 2014) and reflect a plausible range as suggested by global climate models and downscaling experiments in the region (see Echevin et al., 2020 and references therein).

Habitat models

Supervised methods

In the first part of the study, we tuned and compared the two supervised habitat modeling methods, XGBoost and GAM. We also compared different splitting methods (Kfold and shuffle) and evaluated whether the output of the physical-biogeochemical model, CROCO-BioEBUS, can be used as forcing for the habitat models and produce predictions that are comparable with the predictions made by the habitat model when trained with environmental observations.

We fitted a series of XGBoost classifiers to the presence-absence observations with the Python package xgboost (Chen and Guestrin, 2016). To optimize the parameters for XGBoost, we conducted cross-validations together with a brute force grid search for the parameters: learning rate (0.001 to 0.3), depth (2 to 40), and number of estimators (20 to 700). Out of these 200 optimized model configurations, we picked the best model in reference to two metrics: balanced accuracy (in the case of binary classification, this is the area under the ROC curve) and Cohen's Kappa (see Table 2), which are suitable scores for classification problems with imbalanced classes and habitat modeling (Liu et al., 2009; Ben-David, 2008; Wardhani et al., 2019).

We produced Generalized Additive Models (GAM) using the LogisticGAM function, part of the PyGAM package (Table 2), which has a binomial distribution and a logistic link function with the default parameters (Servén and Brummitt, 2018). We adjusted the smoothness level (0.6, 0.1, and 0.0007), number of splines (6, 10, and 20) and the degrees of freedom (1, 2, 3, and 4) using the method described in Section 2.3.2. We then did a cross-validation of the selected model and evaluated train and test metrics for Cohen's Kappa and ROC-AUC (Table 2).

For both XGBoost and GAM models, we used distance from the coast (nm), SST (°C), bathymetry (also referred as water depth; m), depth of the oxycline at 1 ml O 2 L −1 (m) and the month when the sample was taken as predictors. We chose these predictors aiming at maximizing the correlation between observations and output of the CROCO-BioEBUS model while minimizing the correlation between the different predictors (see Figure 3).

The XGBoost and GAM models were cross-validated with 5 splits applying two types of splitting: standard shuffle and Kfold, and we compared the performance of the models (Table 2). Hence, for every step of the cross-validation, 80% of the data was used to train the model and 20% to test. After we had selected the best XGBoost and the best GAM, we applied the models to predict anchovy presence-absence using environmental output from the CROCO-BioEBUS model as predictors instead of observations. We compared the predictions against the observations by calculation ROC-AUC and Kappa again (Table 2). In this way, we could evaluate whether the CROCO-BioEBUS output is close enough to the observations to predict anchovy presence-absence.

Supervised and unsupervised methods comparison

In this section, we compared the supervised methods described above with an unsupervised learning method. For the unsupervised model, we clustered the data cells of the CROCO-BioEBUS model into water types based on their distance from the coast, SST, bathymetry, depth of the oxycline at 1 ml O 2 L −1 and the month when the sample was taken using Gaussian mixture model probability (Figure 4). Then we calculated the probability of finding anchovy in each water type by averaging the presence (1) and absence (0) observation points falling within the water type (Figure 5). Similarly, we calculated the average concentration of anchovies in each water type. Calculating total anchovy biomass is beyond the scope of this paper; therefore, we kept the original units from the raw data of m 2 /nm 2 which is based on acoustic surveys. Since clustering is an unsupervised machine learning method, it is not possible to calculate Kappa and ROC-AUC scores to do cross-validation, so we cannot compare this method to the other models using the standard classification metrics.

In order to compare the three methods, we produced a heatmap of observations averaging at least four data points per grid cell. Then we calculated the mean square root error (RMSE) and correlation coefficient between the grid cells in the observations heatmap and each of our modeled habitat maps. Based on these values, we tuned the parameters of the GAM (see Section 2.3.1) and the number of water types of the clustering (covariance type full) using brute force grid search. A sensitivity of the clustering to the number of water types is available in Appendix 4. All models took the same environmental predictors from CROCO-BioEBUS as input, as described above.

Once we had decided on the best parameters comparing the models with the hindcast observation heatmaps, we applied each of the three models to simulate a climatology of anchovy probability of occurrence taking as predictors the CROCO-BioEBUS output of the sensitivity studies described in Section 2.2.2 to evaluate the climate change impact on their habitat. The maps in our analysis were done using the pcolormesh function in the Matplotlib package (Hunter, 2007). Because the interannual variability in the Northern Humboldt Current System (NHCS) is high in the hindcast simulation that was used to train the models, the comparatively small changes applied to the climatological simulation to represent climate change fall within the range of the training data.

Individual feature contributions to machine-learning predictions were estimated using SHapley Additive exPlanations (SHAP), a game-theoretic framework that decomposes individual predictions into additive feature contributions relative to a base value (Lundberg and Lee, 2017). For the binary classification, SHAP values were interpreted as feature-wise contributions to the probability of the machine learning model, which corresponds in SHAP to log-odds for binary classification. Positive SHAP values increase the anchovy-presence class relative to the base value, whereas negative SHAP values decrease it. The individual classifier output can be reconstructed by adding all feature-specific SHAP values to the base value. In log-odds space, final summed value above zero corresponds to presence, whereas a final summed value below zero corresponds to absence, assuming a default classification threshold of 50% for probabilities. For the spatial interpretation, SHAP values were calculated for each ocean grid cell and subsequently averaged over the respective time period. Global feature relevance was additionally summarized using SHAP beeswarm plots, where each point represents one prediction sample. The horizontal position indicates the magnitude and sign of the feature contribution, while the color represents the corresponding normalized feature value. Thus, the beeswarm plots show both the overall importance of each feature and the distribution of its positive and negative contributions. However, the magnitude of these values are only comparable between features in the same model but not among different models. All calculations were performed separately for 1998 and non-El Niño conditions.

Supervised models training

In our initial approach, we performed the parameter grid search for XGBoost with a standard shuffle splitting method for the maximum depth (3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 30, and 40), number of estimators (20, 40, 50, 70, 100, 150, 200, 250, 300, 350, 400, 450, and 500), and learning rate (0.001, 0.01, 0.1, 0.2, and 0.3), achieving very high train (0.997) and test (0.923) balanced accuracy scores as shown in Table 2. Given the potential for spatial or temporal auto-correlation in dense-sampled observational cruises, these results might be overly optimistic (Salazar et al., 2022). When applying the Kfold splitting method to this complex XGBoost model configuration (maximum depth is 40), the test score dropped by over 30%, strongly indicating model overfitting on likely autocorrelated data. To address this, we also conducted a parameter optimisation with a Kfold splitting method, which splits data into k temporally consecutive folds to prevent overfitting from temporal autocorrelation. Among the models with Kfold-split data, the XGBoost configuration with the highest maximum depth (40) yielded the best training scores. However, the considerably lower test scores suggest overfitting. Conversely, a simpler model with a depth of 3 produced closer test and training scores for both Kappa (0.276 and 0.430, respectively) and the balanced accuracy (0.623 and 0.695, respectively), which were comparable to those of the GAM (Kappa: 0.287, balanced accuracy: 0.629). Therefore, we selected the XGBoost model with maximum depth of 3 as our best model and carried on all further analysis only with this model. In addition, for further validation, we favored Kfold splitting, although the GAM model also showed slightly higher test scores when shuffle-splitted, considering this method is more conservative. Following the integration of data from the CROCO-BioEBUS model instead of observations, the balanced accuracy and Kappa for XGBoost slightly improved from 0.623 to 0.642 and from 0.276 to 0.307, respectively, while for GAM, they decreased from 0.629 to 0.596 and from 0.287 to 0.231, respectively (Table 2). This performance indicates that, despite slight differences, the CROCO-BioEBUS model output is close enough to the observations to produce a reliable habitat map for Peruvian anchovies, demonstrating good generalization across the two different data sets.

Anchovy response to temporal variability

All three models predict a distribution of anchovies near the coast as seen in the observations. In years other than 1998, the distribution extends over the whole Peruvian coastline. In 1998, the probability of finding anchovy decreases with respect to the average of the other years between 1990 and 2010, particularly in the northern part of the domain (Figure 6). In general, the XGBoost model and GAM show a similar map and the same range of probabilities as the occurrence rate in observations. A key difference of these two methods is that, while the XGBoost restrict the regions with the highest probability of occurrence (more than 80%) to the northern and central part of the domain, the GAM model also predicts very high chances of finding anchovies very close to the coast in the southern part of the Peruvian coastline (Figure 6). The clustering map has a narrower range of anchovy probabilities than the XGBoost and GAM, overestimating in regions where the occurrence rate is low and underestimating in regions where finding anchovies is more likely. However, this method exhibited the highest correlation coefficient (0.517) and the lowest RMSE (0.318, Table 2). The correlation coefficient describes the linear relation of the observation and model data, while the RMSE provides the error. On the contrary, the visual comparison suggests that the XGBoost and GAM are better alternatives than the clustering for predicting anchovy probability of occurrence since they show a wider range of predictions reaching zero offshore and up to 1 near the coast. In contrast, clustering represents better the mean state.The clustering method is also able to predict average anchovy concentration in each water type (Figure 5) while XGBoost and GAM were only trained as classifiers. Although concentration and probability of occurrence generally follow a similar pattern, there are differences in their relative changes across water types (Figure 5). The overestimation of the clustering method offshore in contrast to the XGBoost and the GAM is due to the low number of water types. Therefore, offshore water types take part of the signal of anchovy presence in the small areas where they overlap with the high abundance of anchovy near the coast (see Figure 4) although this does not necessarily represent the full extent of the water type.

The spatial distribution of water types classified by the clustering method vary over time (month) and space (Figure 4). Water types 1, 3 and 7 exhibit the highest probability of anchovy occurrence, above 40% (Figures 4, 5), followed by water type 5. These are all coastal water types. Water type 3, which exhibits the highest probability of anchovy occurrence, above 50 %, is found during austral summer months while water type 1, taking the second place corresponds to the winter months (Figures 4, 5). All water types with low probability of anchovy occurrence (below 20%; 0, 2, 6, and 8) extend offshore (Figures 4, 5).

The SHAP analysis shows that depth has a considerably wider range of impacts on the probability of anchovy occurrence in the XGBoost model compared to the other variables, while the clustering and GAM models are more balanced across all the variables (Figure 7). All models generally agree with the direction of the impact from each variable with some slight variations. For example, all models agree that warmer temperatures have a negative impact on anchovies probability of occurrence; however, this relationship is not captured by the clustering method during 1998 (Figure 7).

Anchovy response to climate change

The XGBoost and GAM models predicted a general increase in anchovy probability of occurrence with stronger winds and a decrease in probability with warmer temperatures and weaker winds. A severe maximum decline by about 75% (GAM) and 60% (XGBoost) occurred due to the effect of warmer temperatures (Figure 8). These changes are mostly driven by the increase in temperature in scenario TEMP+ and, to a lesser extent, WIND-, especially in the northern and central parts of the Peruvian coastline. The WIND+ scenario results in a decrease in temperature which correlates with an increase in the probability of anchovy occurrence in the XGBoost and GAM models (Figure 8).

The clustering method showed more ambiguous responses to changes in temperature and wind than the XGBoost and GAM, with parts of the domain where probability of occurrence increases and others where it decreases in the same scenario. The maximum decrease in TEMP+ is less than 15%. Interestingly, there is more variability spatially in the domain than across scenarios. For instance, all scenarios predict a sharp, localized, increase of probability of occurrence at about 10°S. While changes in SST are clearly the main drivers of the response in the XGBoost and GAM models, the clustering model is more strongly affected by the depth of the oxycline (Figure 7). However, this variable shows a mixed response across the domain under the different forcing scenarios, reflecting combined effect of physical and biological responses. In the TEMP+ experiment, increasing temperature strengthens stratification, which generally leads to a shoaling of the oxycline. The WIND+ experiment, though to a lesser degree, also shows a general shoaling pattern. Stronger winds intensify upwelling, lifting isopycnals toward the surface, while the associated increase in biological productivity enhances subsurface oxygen consumption through remineralization, pushing the oxycline shallower still despite the counteracting effect of stronger wind-driven mixing. The WIND- experiment shows the reverse tendency. Compared to SST, which responds relatively quickly and uniformly to the imposed perturbations, the oxycline in the ocean interior generally requires longer adjustment timescales, likely contributing to the patchier spatial pattern of oxycline response across the domain. Overall, the most pronounced changes occur offshore, where anchovies are largely absent. In the southern part of the domain, however, changes in oxycline depth occur closer to the coast, directly affecting the habitat response of anchovies (Figure 8).

The Peruvian anchovy is one of the most studied fisheries due to its economic value but also the complexity of its response to environmental variability. Despite the work done by numerous researchers in the last half of a century, anchovy's behavior is still to some extent unpredictable. In this study, we compared three methods for modeling the habitat of Peruvian anchovies. We used a XGBoost model and a GAM which is a statistical method, both trained via supervised learning. These methods have been commonly used to model habitat. We also tested a new method based on defining water types with unsupervised machine learning, or clustering, and associating anchovy occurrence to each of them. We produced habitat maps with all three models using data from a physical-biogeochemical model and compared their capability for producing realistic maps based on two statistical metrics and a visual assessment of the maps. In Section 4.1 we discuss the implications of using output from physical-biogeochemical models as predictors for habitat models. In Section 4.2 we contextualize the hindcast and climate change predictions from our models. Finally in Section 4.3 we talk about general challenges for modeling the distribution of marine species and how we addressed them.

Applying a physical-biogeochemical model as forcing for habitat models

Producing habitat maps depends on the availability of continuous data over the whole region of interest. This can be obtained, for example, from satellites or in situ observations. However, satellite data is limited to certain variables that can be detected through remote sensing and in situ observations are limited in their spatial coverage. Physical-biogeochemical models, on the other hand, provide a wide range of environmental variables that are not easily observable such as oxycline depth (e.g., Espinoza-Morriberón et al., 2021) and integrated biomass and energy fluxes across the trophic web, which can comprehensively capture the impacts on anchovy feeding and potentially offer a better indicator of fish production (Xue et al., 2024b; Friedland et al., 2012). Furthermore, physical-biogeochemical models can be used to conduct sensitivity experiments (e.g., Echevin et al., 2012; Meier et al., 2019) taking into account the changes in all physical and biogeochemical variables that are relevant for the fish habitat such as temperature and oxygen.

We trained and tested the XGBoost and the GAM models with field observations and then coupled them to a regional physical-biogeochemical model to produce habitat maps. We observed similar values in the balanced accuracy and the Cohen's Kappa score when using predictors from a physical-biogeochemical model in comparison to observational data indicating no loss in predictability due to changing from model data to observations. The XGBoost and GAM models were able to produce very similar maps even when using predictors from a completely different dataset than the one that was used to train them. This indicates that the trained habitat models have high capability to generalize and the relatively low test scores have a different source than the different data domains. Relevant factors could be, for example, the stochastic nature of anchovy appearance (discussed later) or better predictors being needed as model input. However, a similar predictability for both datasets opens the possibility of exploiting the potential of regional physical-biogeochemical models for climate change predictions in the field of habitat modeling.

Earth System Model (ESM) simulations are necessary to predict the impact of climate change on the distribution of marine organisms due to changes in environmental conditions. One way to do this is to add the anomalies simulated by the ESM to the observations used to train the habitat model (e.g., Oliveros-Ramos and Shin, 2023; Sequeira et al., 2014). However, this makes the habitat model susceptible to an unknown error if applied to ESM data since it has been trained and tested only on observations. In addition, the interpolation of global results of the ESM output might not necessarily capture the nuances of a particular region (Echevin et al., 2020). Furthermore, various downscaling and interpolation methods can yield differences in the magnitude of change and in spatial patterns (Pozo Buil et al., 2023). A second method is to use the output of regional models to directly feed habitat models. Nevertheless, the habitat model must be able to generalize well and avoid overfitting to the training data from observations, and be stable enough to handle any noise or bias in the regional model output in order to make good predictions. While regional models also require external forcing from an ESM (e.g., Echevin et al., 2020), some global ESMs such as the Flexible Ocean and Climate Infrastructure (FOCI, Matthes et al., 2020) allow for the implementation of nests with a very fine grid resolution, eliminating the need for downscaling and interpolation at all. Temporal variability in habitat has been suggested to have an impact on how the biomass of marine organisms responds to changes in plankton (Hill Cruz et al., 2022). In addition, fish variability might as well have an impact on the biogeochemistry of the marine ecosystems (Travers-Trolet et al., 2014; Getzlaff and Oschlies, 2017; Hill Cruz et al., 2021; Bianchi et al., 2021). Therefore, being able to run habitat models using direct output from physical-biogeochmeical models opens the possibility for doing dynamic couplings between the two models and capturing the feedbacks between their interaction in end-to-end models. However, before any coupling is attempted, we recommend evaluating which variables in the physical-biogeochemical are close enough to observations to be used as predictors for the habitat model by comparing the performance of the habitat model using in-situ observations and simulated observations (e.g., Tables 1, 2).

Climate change and temporal variability impacts on anchovy habitat

Observations indicate that Peruvian anchovies are generally adapted to cold coastal waters with high productivity and are widespread along the entire Peruvian coast (Bakun and Broad, 2003; Ñiquen Carranza et al., 2000). During El Niño, Peruvian anchovies are found to move closer to the coast and migrate southward, rather than maintaining their typical widespread habitat (Mathisen, 1989; Ñiquen and Bouchon, 2004). This aligns with our findings, where only the very coastal region in the southern part of Peru remains a suitable habitat according to the XGBoost model and GAM during 1998 (Figure 6). The oscillations in the range of the anchovy habitat could potentially influence changes in the fish population; however, as the fish are located further from their typical habitats, they become less accessible for fisheries, leading to an underestimation of fish biomass from commercial catches (Bertrand et al., 2004). Previous studies suggest that environmental changes on different timescales (e.g., interannual and decadal) influence the anchovy habitat in a similar way, which could be indicative of broader climate change impacts (Alheit and Niquen, 2004; Bertrand et al., 2004).

Our clustering method showed that the probability of anchovy occurrence increases in austral summer despite the warmer temperature. In austral winter, anchovies are scattered over a wider area and reach deeper depths (Clark, 1976). This distribution corresponds to our water type 3 (see Figure 4). In contrast, during austral summer, anchovies are found in dense schools near the coast (Clark, 1976), which corresponds to the thin water type 1 in our clustering map (Figure 4). Plankton concentrations near the coast usually increase in summer even though the upwelling is weaker due to a thinner mixed layer depth that decreases light limitation (Echevin et al., 2008; Messié and Chavez, 2015; Xue et al., 2022). Furthermore, Xue and Hill (2025) suggested that a matching in the timing of anchovy hatching and plankton availability further supports the recruitment of their larvae in summer. One caveat of our clustering method is the use of the temporal variable (month) as a continuous variable rather than a factor, generating a gap between December and January that does not exist in reality. However, given that our samples do not cover all months of the year but they are roughly distributed across two seasons (see Appendix Figure 10), and given that we used a relatively small amount of clusters (thirteen), having a continuous month allowed the method to groups months that are close together over the year (with the exception of December to January).

Our climate change experiments indicated a decrease of anchovy occurrence probability along the Peruvian coast with increased temperature and decreased wind in the XGBoost and GAM models (see Figure 8). This is in line with Salvatteci et al. (2022) who suggests a replacement of anchovy by smaller fish as a consequence of climate change. In contrast, increased wind energy may promote habitat expansion for anchovies or counteract the habitat contraction due to warming. Rykaczewski et al. (2015) have suggested that the coastal upwelling zones will be shifted polewards under global warming, resulting in a slightly decreased upwelling strength for the Peruvian coast. However, our simulations did not result in a latitudinal shift of the anchovy habitat for any of the experiments but rather a decrease (TEMP+, WIND-) or increase (WIND+) in the habitat suitability along the entire Peruvian coastline. This is likely a result of our idealized wind speed perturbation, which uses a spatially uniform scaling factor and does not capture potential changes in the wind pattern. Given the uncertainty of projected wind changes in the region, this is a reasonable first assumption.

In this study, the XGBoost and GAM predicted a similar change in both sign and magnitude in the anchovy distribution as a response to climate change. Such methods can be used for moderate climate change projections based on climatological predictors because, in the Northern Humboldt Current System (NHCS), the strong interannual variability observed in the hindcast under present-day climate conditions is stronger than the climate change impacts reflected in the climatology. Hence, we did not attempt to extrapolate in this study. However, our climate change scenarios did not consider the impact of extreme events such as strong El Niño. These events are expected to increase in frequency and intensity in the future (Shin et al., 2022), and this could cause major changes in anchovy habitat and consequently impact fish population and catches (Arias Schreiber et al., 2011).

Despite the general trend observed in models of a preference for colder waters and the implications for climate change, the habitat of Peruvian anchovies remains somewhat unpredictable since they have shown impressive plasticity in recent years. In 2022, the Instituto Público de Investigación de Acuiculture y Pesca (2023) confirmed the presence of anchovies in the Gulf of Guayaquil, in Ecuador, located further north than their usual distribution off Peru. This highlights the adaptable nature of anchovies and the importance of continuous monitoring schemes for ensuring that the most up-to-date changes in species behavior are captured. Since our habitat model was trained with data until 2010, covering only the Peruvian Exclusive Economic Zone, it might not make reliable predictions for the Gulf of Guayaquil. Strong decreases in temperatures were reported in the Easter Pacific Ocean as well as a shallower thermocline in the Gulf of Guayaquil in 2022 during the La Niña event (Instituto Público de Investigación de Acuiculture y Pesca, 2023; Senamhi, 2022), possibly facilitating the occurrence of anchovy further north (Instituto Público de Investigación de Acuiculture y Pesca, 2023). Hence, extending the configuration of our CROCO-BioEBUS model and training the habitat models with up-to-date observations is necessary to be able to simulate an updated hindcast of the anchovy habitat in more recent years and have a clearer picture of their potential response and possibilities of adaptation to climate change. Furthermore, a longer time-series would also allow us to improve the evaluation of both the habitat models developed in this study, as well as the CROCO-BioEBUS model and improve the tuning of key variables that are relevant for the habitat models. The availability of presence-absence data paired with environmental predictors over a long-term is a very valuable tool for improving model performance, especially in face of more recent interannual variability such as El Niño events of 2015–16 and 2017 (Rodŕıguez-Morata et al., 2019). Although the habitat models developed in this project did show a decrease in the probability of occurrence of anchovies during the only strong El Niño event covered in our study period, a longer time-series might further improve their performance under highly changing conditions.

Challenges modeling fish habitat

In order to overcome the challenges of fitting all habitat models used in our study with a common metric and compare both supervised (GAM and XGBoost) and unsupervised (clustering) methods, we opted for producing a heatmap of the observations and comparing the heatmap cells against predicted distribution by the habitat model with input from the physical-biogeochemical model. This method allows us to compare probabilities against occurrence rate rather than single presence-absence points. We were able to train all three models using randomized grid search, to test them against the same set of observations, and to produce maps that provide a realistic distribution for anchovies, concentrating most of them near the coast of Peru and leaving the open ocean free.

A main challenge when modeling marine organisms distributions is the extrapolation of models to regions where observational samples are not available (e.g., Sequeira et al., 2014) or sparse, for instance, the open ocean. This has been handled by adding pseudo-absence data points in areas where it is known that the species of interest does not naturally occur (Luján Paredes, 2016). Our XGBoost model and GAM were trained without pseudo-absence points showing no extrapolation issues, indicating that these models generalize well on unknown data.

A further obstacle when modeling the distribution of highly mobile marine animals such as fish is that their presence in the sampling does not only depend on the suitability of environmental conditions but also on the chances of sampling at the specific time when the school of fish is transiting through the sampled location. Predicting anchovy habitats with a classifier, like the XGBoost model or the GAM, can be challenging if the presence of anchovy schools is probabilistic even under suitable environmental conditions. A key problem is that traditional classifiers might struggle to capture the stochastic nature of anchovy appearances, as they typically aim for deterministic outcomes. This means that, even if the environmental conditions are right, the inherent unpredictability of anchovy locations is not easily learned by standard classifiers, which could lead to inaccuracies in habitat prediction. Furthermore, both XGBoost and GAM require a proper tuning (hyperparamters in the case of XGBoost and splines for GAM). Our clustering method seemed to be robust for handling non-deterministic data and, from a statistical perspective, it was better than the other two methods. It also has potential for predicting abundance which is challenging in datasets that are zero-inflated such as for sparsely distributed organisms (Barry and Welsh, 2002; Lee-Yaw et al., 2022). However, predicting the abundance of anchovy is beyond the scope of this paper so we only evaluated the clustering probability of occurrence predictions. A caveat of the clustering method is that it predicted a very narrow range of probabilities, underestimating in regions where anchovy is common and overestimating in areas where anchovy is absent. If a larger dataset is available, this might be improved by defining more water types. Furthermore, a detailed analysis of its predictions based on every environmental variable showed certain inconsistencies when compared with the supervised methods, which have been extensively used and tested in other studies (e.g., Panzeri et al., 2024; Moore et al., 2009; Hernández-Urcera et al., 2021; Fuster-Alonso et al., 2025; Luan et al., 2021; Grossi et al., 2025; Becker et al., 2020; Effrosynidis et al., 2020; Muhling et al., 2020; Murase et al., 2009). For instance, it gave more importance to the depth of the oxycline while sea surface temperature was a key driver in the other methods. This produced different responses to climate change. Therefore, we recommend caution when using this method. It can support the interpretation of the environmental variables and it produced probability of occurrence predictions that seemed to follow the right pattern in the hindcast simulations. However, further work should be done to understand the behavior of the clustering method and how it uses the environmental variables to make predictions before it can be considered a standard method for habitat modeling.