This study examines the influence of the Indian Ocean Dipole (IOD) on the habitat suitability of Eucheuma cottonii in the coastal waters of Cilacap and Kebumen, which represents the novelty and main objective of this research.
Results indicate that seasonal upwelling in southern Java tends to reduce habitat suitability for E. cottonii.
During normal and positive IOD conditions, stronger upwelling from July to September corresponded with lower Habitat Suitability Index (HSI) values compared with June. In contrast, during the negative IOD phase, when upwelling signals were weaker, HSI values remained relatively stable.
The Kruskal–Wallis test analysis showed significant differences in sea surface temperature (SST), salinity, and nitrate among climate conditions, whereas surface current velocity did not exhibit significant interannual variation.
The combination of sensitivity map analysis and parameter contribution values indicates suggests that bathymetry (water depth) may play a relatively important role in shaping habitat suitability patterns, while other oceanographic parameters (SST, salinity, nitrate, and surface current velocity) function as dynamic variables whose influence may vary depending on the prevailing IOD climate conditions.
E. cottonii is one of the major seaweed commodities cultivated in Indonesia (Manurung et al., 2021; Wahyuni et al., 2023; Wijayanto et al., 2020). The economic benefits obtained from its cultivation can reach IDR 19,500,000 (US$1184.59) per cultivated area or IDR 23,900,000 (US$1451.88) per hectare per harvest season (Wahyuni et al., 2023), with a cultivation period ranging from 30 to 65 days (Wijayanto et al., 2020). The success of seaweed cultivation is strongly influenced by site selection (Casadebaig et al., 2022) therefore, identifying suitable cultivation areas is crucial. One rapidly developing approach for this purpose is habitat suitability modeling (Bertelli et al., 2022; Stephenson et al., 2021; Stuart et al., 2021). Habitat suitability in coastal waters is largely affected by climate variability and oceanographic factors (Gokturk et al., 2022; Sun et al., 2024).
The coastal upwelling system along the southern coast of Java plays a significant role in determining marine productivity and controlling oceanographic conditions in the region (Napitupulu, 2025; Rochmatika & Bahtiar, 2023). In general, upwelling is defined as the upward movement of subsurface waters to the surface, driven by Ekman pumping resulting from persistent wind forcing and the Coriolis effect (Koropitan et al., 2021). In addition to the influence of the southeast monsoon, coastal upwelling along southern Java is frequently modulated by large-scale climate phenomena such as the El Niño–southern oscillation (ENSO) and the Indian Ocean Dipole (IOD) (Widagdo et al., 2025), which can either intensify or weaken the upwelling system (Hafiz et al., 2024; Oktaviani et al., 2021). Consequently, the coastal waters of Cilacap and Kebumen, as part of the southern Java upwelling system, are directly affected by these climate phenomena.
Coastal upwelling along southern Java typically occurs during the peak of the southeast monsoon between June and September (JJAS) (Budiman et al., 2022; Horii et al., 2023; Koropitan et al., 2021). The IOD index data obtained from the Bureau of Meteorology (BOM) indicate that positive IOD conditions occurred during JJAS in 2019, negative IOD conditions occurred during JJAS in 2016, while 2018 can be classified as a normal year during the same period.
The IOD is a climate phenomenon characterized by anomalous cooling or warming of sea surface temperature (SST) in the eastern or western Indian Ocean. A positive IOD event occurs when SST in the eastern Indian Ocean (western Sumatra and Java) cools while SST in the western Indian Ocean (eastern Africa) warms, whereas the opposite pattern characterizes negative IOD events (Koropitan et al., 2021; Wang et al., 2024). This near-consecutive interannual climate variability (positive IOD in 2019, negative IOD in 2016, and normal conditions in 2018) provides a valuable opportunity to investigate the habitat suitability of E. cottonii in relation to upwelling dynamics and IOD variability in the coastal waters of Cilacap and Kebumen, Indonesia.
One approach that can be applied to address this issue is habitat suitability modeling. Several studies have linked habitat suitability modeling with IOD variability. Lan et al. (2013) examined the habitat suitability of yellowfin tuna (Thunnus albacares) in the western Indian Ocean and reported a significant increase in habitat suitability during negative IOD events and a decrease during positive IOD events. Lan et al. (2015) investigated the habitat suitability of swordfish (Xiphias gladius) in the Indian Ocean, showing significant increases along the northwestern Indian Ocean coast during positive IOD events and in the southern Indian Ocean during negative IOD events. Koropitan et al. (2021) analyzed the habitat suitability of mackerel tuna (Euthynnus affinis) influenced by negative IOD conditions along the southern Java coast and found a significant decline in habitat suitability during negative IOD events. Wen et al. (2025) studied the habitat suitability of squid (Sthenoteuthis oualaniensis) in the northern Indian Ocean and demonstrated that suitable habitat areas were smaller during positive IOD events and larger during negative IOD events. To date, no studies have explicitly modeled seaweed aquaculture habitat suitability while simultaneously examining the influence of the IOD.
A review of the literature over the past decade indicates that most habitat suitability studies of E. cottonii remain largely static, typically conducted for a single time period and without integrating physical oceanographic dynamics or regional climate variability. Examples include studies conducted in Mandar Bay, West Sulawesi by Rusdi (2018); Tidung Island by Utama and Handayani (2018); the eastern coast of Tarakan Island by Lestari et al. (2019); Pasiea District, North Buton by Salihin et al. (2019); Parepare Bay by Damis and Saenong (2020); Ambon Baguala Bay by Lase et al. (2020); Sarawandori Bay, Yapen Papua by Numberi et al. (2020); West Sorkam, Central Tapanuli by Manurung et al. (2021); Takalar Lama by St Madina et al. (2022); Majene waters by Arbit et al. (2024); and Lapang Island waters by Wabang and Plaimo (2024). Most of these studies primarily focused on mapping habitat suitability based on local environmental parameters without linking them to larger-scale oceanographic forcing mechanisms.
To date, only Lestari et al. (2019) have related E. cottonii habitat suitability to climate variability associated with the ENSO, demonstrating that the potential cultivation area was larger during La Niña compared to El Niño conditions. However, the influence of major oceanographic processes such as coastal upwelling and the role of variability associated with the IOD on E. cottonii habitat suitability has not yet been explicitly investigated, particularly in the southern coastal waters of Java, a region known for its strong seasonal upwelling dynamics.
The main novelty of this study lies in integrating upwelling processes and IOD variability into the habitat suitability analysis of E. cottonii. This approach not only describes the spatial patterns of habitat suitability but also explains the oceanographic and climatic mechanisms that control its variability. This study represents the first attempt to explicitly assess the combined influence of upwelling and the IOD on E. cottonii habitat suitability in the coastal waters of Cilacap and Kebumen. The findings are expected to provide a stronger scientific basis for climate-adaptive seaweed aquaculture management and contribute to advancing theoretical understanding of how IOD-related climate variability influences coastal aquaculture environments.
This study was conducted in the coastal waters of Cilacap and Kebumen, Indonesia, located between 108.5°E–109.8°E and 7.3°S–8.0°S (Figure 1).
Study area map.
Environmental parameters included SST, salinity, surface current velocity (zonal [u] and meridional [v]), and nitrate concentration. Oceanographic data were obtained from the Copernicus Marine Service reanalysis products. SST, salinity, and surface current velocity were derived from the product GLOBAL_MULTIYEAR_PHY_001_030, which has a spatial resolution of 0.083° × 0.083°. Nitrate concentration data were obtained from the product GLOBAL_ANALYSISFORECAST_BGC_001_028, which has a spatial resolution of 0.25° × 0.25°. Both datasets provide monthly averaged fields that have undergone operational quality control and data assimilation procedures. Bathymetric data were obtained from the General Bathymetric Chart of the Oceans (GEBCO) dataset, specifically the GEBCO 2025 Grid. The bathymetric dataset has a spatial resolution of 0.0041667° × 0.0041667° (equivalent to 15” or approximately 450 m at the equator). The analysis focused on the JJAS period to represent three contrasting climate conditions: negative IOD (2016), normal conditions (2018), and positive IOD (2019). Seasonal means were calculated as the arithmetic average of the monthly values during the JJAS period. All datasets were processed in NetCDF format.
Oceanographic variables from Copernicus Marine Service are operationally validated through the integration of numerical models, satellite observations, and in situ measurements. Given the quality control procedures implemented by the data provider, additional independent validation was not performed in this study. Nevertheless, the interpretation of results considers the spatial resolution limitations of the datasets (approximately 0.083°–0.25°), particularly in complex coastal environments. The GEBCO bathymetry grid represents a global compilation derived from hydrographic surveys and satellite altimetry and has undergone international quality control procedures. In this study, bathymetric data were primarily used as a depth constraint for habitat suitability filtering rather than for detailed seafloor morphology analysis; therefore, additional field validation was not conducted. However, it is acknowledged that local bathymetric accuracy may vary depending on the density of surveys in specific regions.
A habitat suitability analysis model for E. cottonii aquaculture was developed using a Boolean threshold approach within a geographic information system (GIS) framework and implemented using the Python programming language. The Boolean method assigns a score of 1 to environmental parameters that fall within suitable ranges and a score of 0 to those that fall outside suitable ranges (Ukhti et al., 2021). The suitability thresholds for E. cottonii cultivation, based on Manurung et al. (2021), were defined as follows: SST of 26°C–33°C, salinity of 25–35 PSU, nitrate concentration of 0.1–4.4 mg L−1 (equivalent to 1.61–70.96 mmol m−3), current velocity of 0.1–0.4 m s−1, and bathymetry (water depth) of 0.3–10 m. After scoring each parameter, the Habitat Suitability Index (HSI) was calculated using the Arithmetic Mean Model (AMM) (Wen et al., 2025), expressed as:
HSI = Habitat Suitability Index.
Si = suitability score for the ith parameter.
n = number of parameters.
This approach does not assume that all environmental parameters exert identical biological effects. Instead, it aims to avoid assigning subjective parameter weights when quantitative information regarding the relative contribution of each variable is not available. Previous studies indicate that the most influential environmental parameters affecting the growth of E. cottonii vary across locations. For instance, Manurung et al. (2021) identified current velocity as the dominant factor, whereas Numberi et al. (2020) emphasized water transparency as the primary parameter, and Lestari et al. (2019) highlighted the importance of dissolved oxygen (DO). Other studies suggest that multiple parameters may act simultaneously. Utama and Handayani (2018) reported substrate type, current velocity, and salinity as key factors; Wabang and Plaimo (2024) emphasized substrate, wave exposure, currents, water clarity, and depth; Arbit et al. (2024) identified temperature, salinity, and nutrients (nitrate and phosphate) as influential variables; while Salihin et al. (2019) underlined the role of salinity, transparency, phosphate, and nitrate. The variability of these findings indicates that no single environmental parameter consistently dominates across regions. Consequently, assigning differential weights without local biological evidence may introduce subjective bias. Therefore, equal weighting through a Boolean approach was adopted as a conservative and transparent strategy, allowing all ecologically relevant parameters to contribute equally to the habitat suitability assessment.
The compensatory assumption embedded in this model reflects the ecological tolerance range of E. cottonii, whereby suboptimal conditions in one parameter may still be compensated by favorable conditions in others, provided that ecological thresholds are not exceeded. This approach is appropriate for studies aiming to explore spatial patterns of habitat suitability rather than to quantitatively predict aquaculture production. Under conditions of limited aquaculture production data, the Boolean approach combined with arithmetic averaging provides a reproducible method with minimal subjective assumptions and is supported by the empirical variability reported in previous studies.
All analyses were conducted using the Python programming language with the primary libraries xarray, numpy, scipy, matplotlib, and cartopy. The Python scripts used in this study are provided as Supplementary Material to ensure analytical reproducibility.
The pre-processing stage began by loading each dataset using xarray.open_dataset() function and extracting the relevant variables. For variables containing a depth dimension, only the surface layer (depth = 0 m) was selected using the isel (depth = 0) function, as the study focuses on habitat suitability for seaweed cultivation in the upper water column. All variables were subsequently aligned to the spatial grid and temporal dimension of SST as the reference dataset using the interp_like() interpolation method with a nearest-neighbor approach. This procedure ensures spatial and temporal consistency among all environmental parameters prior to further analysis. Bathymetric data were also interpolated to the SST grid using the same method to ensure identical coordinate systems and spatial resolution across all layers. Surface current velocity was calculated from the zonal and meridional current components using the equation:
u = represents the zonal current component and
v = represents the meridional current component.
A threshold-based suitability selection was then applied to each environmental parameter using ecological thresholds derived from previous literature: SST (26°C–33°C), salinity (25–35 PSU), nitrate (1.61–70.96 mmol m−3), current velocity (0.1–0.4 m s−1), and water depth (0.3–10 m). Each parameter was converted into a binary suitability score (0 = unsuitable; 1 = suitable). The final suitability score was calculated as the arithmetic mean of the five parameters, resulting in a HSI ranging from 0 to 1. To reduce spatial noise caused by grid discretization, Gaussian smoothing was applied using the gaussian_filter function with a sigma value of 1. The resulting suitability map was then spatially cropped to the study area covering the coastal waters of Cilacap and Kebumen (108.5°E–109.8°E; 8.0°S –7.3°S) using coordinate-based slicing. The analyzed coastal grid cells implicitly represent marine pixels that satisfy the depth criterion and contain valid values for all environmental parameters after interpolation and masking procedures. Spatial visualization was produced using the Plate Carrée projection through the cartopy library, and all outputs were exported as raster images (.png) for each analyzed month.
Although the physical datasets from the Copernicus Marine Service and bathymetric data from the GEBCO have undergone operational validation and therefore do not require additional independent validation, the HSI model developed in this study has not been validated against in situ observations or actual seaweed aquaculture production data. Therefore, the results represent potential habitat suitability based on biophysical conditions rather than empirical validation of cultivation performance in the field. This limitation is explicitly acknowledged and highlights the need for future studies integrating field survey data or aquaculture production records to improve the predictive accuracy of habitat suitability models.
Differences in environmental parameters (SST, salinity, nitrate concentration, and surface current velocity) among climate conditions (negative IOD in 2016, neutral conditions in 2018, and positive IOD in 2019) during the JJAS period were evaluated using a non-parametric statistical approach. The analyzed data represent the midpoint values of each parameter within the study area for each month during JJAS, such that each climate condition is represented by four monthly observations. Initial visualization was performed using boxplots to illustrate the median, interquartile range (IQR), and variability among climate groups. Given the limited sample size, individual data points were overlaid on the boxplots using scatter markers to enhance transparency of the underlying data distribution. All visualizations were generated using the Matplotlib library in Python.
Hypothesis testing was conducted using the Kruskal–Wallis test, implemented through the kruskal function in the SciPy statistical module. This test was selected because of the relatively small sample size and its ability to compare multiple independent groups without assuming normal data distribution or homogeneity of variance. The Kruskal–Wallis test is a rank-based method that evaluates whether the medians of more than two independent groups differ significantly. The test statistic is expressed as the H value, with statistical significance evaluated at α = 0.05. A p-value below 0.05 was interpreted as evidence of a statistically significant difference in parameter distributions among climate conditions. Because the analysis is based on spatially aggregated midpoint values and a limited temporal sample, statistical outcomes were interpreted cautiously as exploratory comparative evidence rather than causal inference. All analyses were performed using Python, and the complete scripts are provided as supplementary material to ensure transparency and reproducibility of the study.
Sensitivity analysis was conducted to evaluate the influence of each environmental parameter on the developed HSI model. This analysis is important in habitat modeling because the results of habitat suitability models can be influenced by the quality of input data, the selection of environmental variables, and the configuration of the model used. Therefore, sensitivity evaluation is necessary to ensure that the model produces stable and reliable habitat estimates (Arsenault et al., 2025). In this study, the sensitivity analysis was performed using a leave-one-parameter-out approach, also known as jackknife sensitivity analysis. This approach is widely used in species distribution and habitat modeling studies to identify environmental variables that have the strongest influence on model outcomes. The principle of this method is to recalculate the model output after removing one environmental parameter and then compare the resulting model output with the original model that includes all parameters (Wan et al., 2024).
The first step was to calculate the full HSI using all environmental parameters included in the model, namely SST, salinity, nitrate, current velocity, and bathymetry (water depth). The model was then recalculated by removing one environmental parameter at each iteration, producing an HSI value without the i-th parameter. The change in the habitat suitability index resulting from the removal of each parameter was calculated as the sensitivity value using the following equation:
HSIfull = the habitat suitability index calculated using all environmental parameters
HSI(−i) = the habitat suitability index calculated after the i- th parameter is removed from the model.
This sensitivity value represents the magnitude of change in the model when a particular parameter is excluded from the HSI calculation. A larger change indicates a stronger influence of that parameter on the habitat suitability model. In addition, the relative contribution of each parameter to the model was calculated using the following equation:
The relative contribution value is used to identify the environmental parameters that are most dominant in determining habitat suitability. Parameters with higher contribution values are considered to have a greater influence on the variability of HSI values.
This sensitivity analysis was conducted under three different climate conditions, namely normal conditions, negative IOD, and positive IOD, allowing an evaluation of how the contribution of environmental parameters to the habitat suitability model varies under different climate variability scenarios. The Python scripts used for the sensitivity analysis are provided in the supplementary materials.
The spatial distribution of SST during the JJAS period for 2016 (negative IOD), 2018 (normal conditions), and 2019 (positive IOD) along the southern coast of Java is shown in Figure 2. The results indicate that the IOD phenomenon strongly influences the spatial extent and intensity of coastal upwelling along southern Java. Upwelling is characterized by the presence of low SST values, reflecting the upward movement of cooler subsurface waters to the surface (Koropitan et al., 2021; Rachman et al., 2024). Under normal conditions in 2018, upwelling began to develop along the southern coast of Java in July, as indicated by the emergence of SST values below 25°C (blue shading). The upwelling persisted until September, reaching its peak intensity in August 2018. During positive IOD conditions in 2019, upwelling was observed from JJAS, with a pronounced peak in August and a noticeably larger spatial extent compared to normal conditions in 2018. In contrast, during negative IOD conditions in 2016, upwelling was not detected throughout the JJAS period. Instead, SST values were significantly higher than those observed during the normal conditions in 2018. These findings are consistent with previous studies, which reported that positive IOD events enhance upwelling intensity, whereas negative IOD events suppress upwelling along the southern coast of Java (Rachman et al., 2024).
Spatial distribution of sea surface temperature (SST) along the southern coast of Java during June, July, August, and September under negative IOD conditions (2016), normal conditions (2018), and positive IOD conditions (2019).
The results of the habitat suitability model for E. cottonii aquaculture are presented in Figure 3. The model produces HSI values ranging from 0 to 1, where higher values indicate more suitable environmental conditions for E. cottonii cultivation. The results demonstrate that the occurrence of coastal upwelling exerts a negative influence on the habitat suitability model for E. cottonii aquaculture. Under normal conditions in 2018, the development of upwelling during July to September resulted in a decrease in HSI values, from a dominant value of 0.4 in June–0.2 in July and approximately 0.3 in August and September. Although upwelling is known to enhance primary productivity in coastal waters (Rachman et al., 2024), this increase was not accompanied by an improvement in HSI values for E. cottonii aquaculture.
Habitat suitability model for Eucheuma cottonii aquaculture in the coastal waters of Cilacap and Kebumen under negative IOD (2016), normal (2018), and positive IOD (2019) conditions during June, July, August, and September.
A similar pattern was observed during positive IOD conditions in 2019, when increased upwelling intensity did not lead to improved habitat suitability. During this year, HSI values were predominantly 0.4 in June and declined to approximately 0.3 from July to September. In contrast, during negative IOD conditions in 2016, when upwelling was absent in the coastal waters of Cilacap and Kebumen, the habitat suitability model for E. cottonii aquaculture remained stable, with HSI values predominantly around 0.4 throughout the JJAS period (Figure 3). Based on Figure 3, the most suitable areas for E. cottonii cultivation, characterized by the highest HSI values across all months, are located in the coastal waters of Cilacap, particularly in the Segara Anakan region and along Teluk Penyu Beach (108.8°E–109.05°E). Additional suitable areas are identified in the eastern part of Kebumen waters, specifically from Suwuk Beach to Laguna Lembupuro Beach, within the longitude range of 109.42°E–109.8°E.
The HSI values obtained in this study ranged from 0 to 0.4 during the observation period. It should be emphasized that the model applied in this study follows a Boolean threshold approach based on five environmental parameters. Consequently, the HSI value represents the proportion of parameters that simultaneously meet the defined ecological suitability thresholds. Under this formulation, an HSI value of 0.4 indicates that two out of the five environmental parameters fall within their optimal ranges at a given grid cell. Within the framework of the classical Habitat Suitability Index developed under the Habitat Evaluation Procedures, the maximum value (HSI = 1) is only achieved when all environmental variables simultaneously reach their optimal conditions. Therefore, the absence of values ≥0.6 in this study indicates that no location satisfied at least three optimal environmental parameters simultaneously during the JJAS period (Zhang et al., 2025).
Habitat suitability classification in this study does not adopt universal thresholds (e.g., ≥0.6 as ‘high suitability’), because habitat modeling literature highlights that classification boundaries are context-dependent and strongly influenced by model structure and analytical objectives (Kroth et al., 2025). Given the binary threshold approach and the regional-scale environmental datasets used in this study, HSI values are interpreted relatively, focusing on spatial and temporal differences among climate conditions rather than defining absolute suitability categories (Suresh et al., 2025). Accordingly, an HSI value of 0.4 should be interpreted as a relatively moderate level of suitability within the comparative framework of this study, rather than as an indication that the habitat is biologically optimal.
The upwelling process, which transports subsurface waters to the surface, not only alters SST but also modifies surface salinity, current velocity, and nutrient concentrations (Awo et al., 2022; Budiman et al., 2022; Li et al., 2024). The SST range used to construct the habitat suitability model for E. cottonii aquaculture was 26°C–33°C (Manurung et al., 2021). Under normal conditions in 2018, SST in the coastal waters of Cilacap and Kebumen ranged from 27.3°C to 27.7°C in June, followed by a decrease in July to 24.9°C–25.4°C, 24.4°C–25.1°C in August, and 24.9°C–25.7°C in September (Figure 4). This decline in SST was associated with the occurrence of upwelling, resulting in SST values falling below the suitable range for E. cottonii cultivation during the upwelling period. During positive IOD conditions in 2019, enhanced upwelling intensity further reduced SST values, causing SST to deviate more strongly from the suitable range during July to September. SST values ranged from 25.9°C to 26.6°C in June, 24.7°C–25.3°C in July, 23.7°C–24.3°C in August, and 23.7°C–24.8°C in September. These results are consistent with Iskandar et al. (2022), who reported a significant SST decrease along the southern coast of Java due to upwelling processes during June–October 2019.
Spatial distribution of sea surface temperature (SST) in the coastal waters of Cilacap and Kebumen under negative IOD (2016), normal (2018), and positive IOD (2019) conditions during June, July, August, and September.
In contrast, negative IOD conditions in 2016 were characterized by significantly higher SST values compared to normal conditions, with SST ranging from 29.8°C to 30.3°C in June, 29.4°C–29.8°C in July, 29.2°C–29.6°C in August, and 29.6°C–29.8°C in September. These conditions resulted in SST values that remained within the suitable range for E. cottonii cultivation throughout the JJAS period. Koropitan et al. (2021) reported that during July–September 2016, monthly SST exceeded 29°C across parts of the western Sumatra coast and the southern coast of Java. This SST increase is attributed to negative IOD conditions, which extended their influence into the Java Sea and the southern Java upwelling system.
The salinity range used in the habitat suitability model for E. cottonii aquaculture was 25–35 PSU (Manurung et al., 2021). The salinity distribution during the JJAS period under negative IOD conditions (2016), normal conditions (2018), and positive IOD conditions (2019) remained entirely within the suitable range for seaweed cultivation. Under normal conditions in 2018, salinity values during JJAS ranged from 34.1 to 34.3 PSU, 34.5–34.6 PSU, 34.7–34.8 PSU, and 34.6–34.8 PSU, respectively. During negative IOD conditions in 2016, salinity values ranged from 33.2 to 33.4 PSU, 33.5–33.7 PSU, 33.7–33.8 PSU, and 33.3–34.5 PSU. During positive IOD conditions in 2019, salinity values ranged from 34.2 to 34.6 PSU, 34.6–34.7 PSU, 34.7–34.8 PSU, and 34.6–34.7 PSU (Figure 5).
Spatial distribution of surface salinity in the coastal waters of Cilacap and Kebumen under negative IOD (2016), normal (2018), and positive IOD (2019) conditions during June, July, August, and September.
The increase in salinity associated with upwelling was relatively small, at approximately 0.2 PSU, while the occurrence of positive IOD conditions resulted in an additional increase of only about 0.1 PSU compared to normal conditions. In contrast, negative IOD conditions were associated with a decrease in surface salinity of up to 1 PSU. Huang et al. (2025) reported that during positive IOD events, equatorial easterly winds induce upwelling, strengthening anticyclonic salt advection, and leading to positive salinity anomalies. Chen et al. (2022) reported that during negative IOD events, surface salinity in the Indian Ocean can decrease by up to 2 PSU due to enhanced convergence, which increases freshwater transport.
Although salinity did not exhibit sufficient variation to fall outside the optimal range (25–35 PSU) throughout the entire observation period, equal weighting was retained in this model based on conceptual and methodological considerations. The HSI model employed was designed to represent suitability conditions simultaneously rather than to measure the relative sensitivity of each parameter to specific temporal variations. Thus, the stability of salinity within the optimal range does not imply that this variable is ecologically unimportant; rather, it indicates that during the JJAS period in the study area, salinity does not constitute a primary limiting factor for the cultivation of E. cottonii. Physiologically, tropical seaweeds are known to possess a relatively narrow salinity tolerance, and even minor changes can affect osmoregulatory processes and growth (Aris et al., 2021; Tilaar et al., 2025). Therefore, although interannual variation ranged only from 0.1 to 1 PSU (Figure 6), this parameter was retained to ensure that the model remains sensitive to potential extreme conditions beyond the analyzed period. Furthermore, habitat modeling literature emphasizes that assigning differential weights without a robust local quantitative basis may introduce subjective bias (Arsenault et al., 2025).
Boxplot distribution of environmental parameter under different climate conditions.
The nitrate concentration range used in the habitat suitability model for E. cottonii aquaculture was 0.1–4.4 mg L−1, equivalent to 1.61–70.96 mmol m−3. Under normal conditions in 2018, nitrate concentrations in the coastal waters of Cilacap and Kebumen were relatively low in June, at approximately 0.4 mmol m−3. The onset of upwelling during July, August, and September increased nitrate concentrations to approximately 4.0, 4.6–5.0, and 3.8 mmol m−3, respectively (Figure 7). Rachman et al. (2024) reported that upwelling processes transport nutrient-rich subsurface waters to the surface, thereby enhancing primary productivity in coastal waters.
Spatial distribution of nitrate concentrations in the coastal waters of Cilacap and Kebumen under negative IOD (2016), normal (2018), and positive IOD (2019) conditions during June, July, August, and September.
During positive IOD conditions in 2019, nitrate concentrations increased markedly, reaching 3.5–4.0 mmol m−3 in June, approximately 4.5 mmol m−3 in July, 8.1 mmol m−3 in August, and exceeding 10 mmol m−3 in September. In contrast, during negative IOD conditions in 2016, nitrate concentrations decreased substantially, reaching approximately 0.2 mmol m−3 in June, 0.1 mmol m−3 in July, 0 mmol m−3 in August, and 0.1 mmol m−3 in September. Maysarah et al. (2022) demonstrated that positive IOD events are significantly associated with increased nitrate concentrations, particularly during the June to December period in southern Java waters, whereas the opposite pattern occurs during negative IOD events.
In the habitat suitability model applied in this study, a threshold-based Boolean approach was used, whereby any environmental parameter falling outside the defined optimal range, either below the lower limit or exceeding the upper limit, was automatically assigned a score of 0 (unsuitable). Therefore, during positive IOD conditions, when nitrate concentrations exceeded the upper threshold of the optimal range (Manurung et al., 2021), these values were classified as unsuitable for the nitrate parameter. Although nitrate enrichment associated with upwelling can enhance primary productivity, excessively high nutrient concentrations do not necessarily benefit the growth of cultivated macroalgae, as they may induce physiological imbalance and intensify competition with phytoplankton and epiphytic organisms (Hurd et al., 2014). Accordingly, in the context of E. cottonii cultivation, nitrate levels exceeding the optimal range were treated as suboptimal conditions within the habitat suitability assessment framework used in this study.
The current velocity range used in the habitat suitability model for E. cottonii aquaculture was 0.1–0.4 m s−1 (Manurung et al., 2021). Under normal conditions in 2018, current velocities in the coastal waters of Cilacap and Kebumen ranged from 0.04 to 0.20 m s−1 in June, 0.07–0.21 m s−1 in July, 0.03–0.23 m s−1 in August, and 0.26–0.36 m s−1 in September. The occurrence of IOD events, whether positive or negative, did not substantially alter current velocities in the coastal waters of Cilacap and Kebumen. During positive IOD conditions in 2019, current velocities ranged from 0.02 to 0.05 m s−1 in June, 0.12–0.15 m s−1 in July, 0.12–0.21 m s−1 in August, and 0.13–0.18 m s−1 in September. Similarly, during negative IOD conditions in 2016, current velocities ranged from 0.03 to 0.12 m s−1 in June and July, 0.10–0.18 m s−1 in August, and 0.03–0.10 m s−1 in September (Figure 8).
Spatial distribution of surface current velocity in the coastal waters of Cilacap and Kebumen under negative IOD (2016), normal (2018), and positive IOD (2019) conditions during June, July, August, and September.
These results differ from previous studies. Wijaya et al. (2024) reported that the south Java current (SJC) exhibits a negative relationship with IOD and ENSO during the June–August period, while Maysarah et al. (2022) found that the SJC strengthens during negative IOD events and weakens during positive IOD events. The discrepancy between the present results and previous findings is likely attributable to the relatively small spatial extent of the study area and its proximity to the coastline (≤8°S). Under these conditions, the influence of the SJC and large-scale climate variability on observed current velocities may be diminished, resulting in weaker detectable signals in the reanalysis data.
The bathymetry (water depth) range used in the habitat suitability model for E. cottonii aquaculture was 0.3–10 m (Manurung et al., 2021). As shown in Figure 9, suitable depth conditions were primarily confined to the eastern part of Kebumen waters, particularly within the longitude range of 109.35°E–109.8°E. Limited suitable areas were also identified in the coastal waters of Cilacap, specifically around the Segara Anakan region and Teluk Penyu Beach (Figure 9). These areas represent shallow coastal zones that meet the depth requirements for E. cottonii cultivation and largely determine the spatial distribution of high habitat suitability values observed in the habitat suitability model.
Bathymetric distribution in the coastal waters of Cilacap and Kebumen.
Among all environmental parameters used to construct the habitat suitability model, the spatial and temporal variation of SST (Figure 4) exhibits the strongest similarity to the spatial distribution and temporal variation of the HSI (Figure 3). This finding indicates that SST plays a dominant role in controlling the dynamics of habitat suitability for E. cottonii aquaculture.
The occurrence of upwelling, which lowers SST in the coastal waters of Cilacap and Kebumen, was found to negatively affect the habitat suitability of E. cottonii (Figure 3). This response contrasts with that of many other marine resources, whose productivity typically increases during upwelling events. These results raise an important question regarding the optimal timing of E. cottonii cultivation based on habitat suitability conditions, which should be addressed in future studies to support more effective cultivation recommendations.
Positive IOD events, which are associated with SST cooling in the coastal waters of Cilacap and Kebumen, did not exert a significant influence on the habitat suitability of E. cottonii. In contrast, negative IOD events, characterized by SST warming, had a significant positive effect on habitat suitability compared to normal conditions. This result highlights the critical sensitivity of E. cottonii cultivation to thermal conditions and suggests that warmer SSTs associated with negative IOD events provide more favorable environmental conditions for seaweed aquaculture in the study area. These results are consistent with previous findings by Lestari et al. (2019), who reported that during La Niña conditions characterized by relatively warmer SST and weaker upwelling, the potential area suitable for seaweed cultivation tends to be larger than during El Niño periods, which are typically associated with intensified upwelling.
Monthly midpoint values were used to represent environmental conditions and were treated as independent observations in the Kruskal–Wallis test to compare differences among climate conditions (negative IOD, normal, and positive IOD). The ranges and midpoint values of monthly environmental parameters under different climate conditions in the coastal waters of Cilacap and Kebumen are presented in Table 1, while the distribution of each parameter is illustrated using boxplots in Figure 6. The distribution of environmental parameters shows distinct responses to climate conditions. SST exhibits a clear decrease during positive IOD conditions, while nitrate concentration increases substantially, reflecting the intensification of coastal upwelling processes. In contrast, salinity and surface current velocity show relatively small variations among climate conditions, suggesting that these parameters are more strongly influenced by local factors than by regional climate variability. Because the environmental parameter data were not assumed to follow a normal distribution, a non-parametric approach using the Kruskal–Wallis test was considered the most appropriate method for evaluating differences among climate conditions.
Range and midpoint values of monthly environmental parameters under different climate conditions in the coastal waters of Cilacap and Kebumen.
| Year/Condition | Month | Value range | Mean value |
|---|---|---|---|
| SST (°C) | |||
| Negative IOD (2016) | June | 29.8–30.3 | 30.05 |
| Normal (2018) | June | 27.3–27.7 | 27.50 |
| Positive IOD (2019) | June | 25.9–26.6 | 26.25 |
| Salinity (PSU) | |||
| Negative IOD (2016) | June | 33.2–33.4 | 33.30 |
| Normal (2018) | June | 34.1–34.3 | 34.20 |
| Positive IOD (2019) | June | 34.2–34.6 | 34.40 |
| Nitrate (mmol m−3) | |||
| Negative IOD (2016) | June | 0.2 | 0.20 |
| Normal (2018) | June | 0.4 | 0.40 |
| Positive IOD (2019) | June | 3.5–4.0 | 3.75 |
| Current velocity (m s−1) | |||
| Negative IOD (2016) | June | 0.03–0.12 | 0.075 |
| Normal (2018) | June | 0.04–0.20 | 0.120 |
| Positive IOD (2019) | June | 0.02–0.05 | 0.035 |
IOD, Indian Ocean Dipole; SST, sea surface temperature.
The Kruskal–Wallis test results indicate significant differences in SST (p = 0.018 <0.05), salinity (p = 0.024 <0.05), and nitrate concentration (p = 0.018 <0.05) among climate conditions. In contrast, surface current velocity does not show a significant difference (p = 0.524 >0.05) (Table 2). In the Kruskal–Wallis test, the H statistic represents the magnitude of distributional differences among groups in this case, climate conditions (negative IOD, normal, and positive IOD) (Wilcox, 2012). Larger H values accompanied by significant p-values indicate a stronger influence of the differentiating factor, namely climate conditions (Rustin & Arifin, 2026; Saensouk et al., 2025). In this study, the H values obtained from the Kruskal–Wallis test are as follows: SST (H = 8.00), salinity (H = 7.45), nitrate concentration (H = 8.03), and surface current velocity (H = 1.29) (Table 2). The significant differences in SST, salinity, and nitrate concentration among climate conditions reflect the strong influence of basin-scale climate variability, particularly the IOD, on the oceanographic dynamics of the Cilacap and Kebumen coastal waters. During positive IOD conditions, the significant cooling of SST indicates an intensification of coastal upwelling along the southern coast of Java. This process promotes the upward transport of colder, more saline, and nutrient-rich subsurface water toward the surface, which directly explains the observed increase in nitrate concentrations during this period. Conversely, the opposite tendency occurs during negative IOD conditions (Awo et al., 2022; Budiman et al., 2022; Li et al., 2024). The absence of significant differences in surface current velocity among climate conditions suggests that surface current variability in this region is primarily controlled by local and regional factors, such as coastal topography, cross-shore circulation patterns, and intraseasonal wind fluctuations, rather than by interannual climate signals (Gu & Mao, 2024). This indicates that coastal environmental responses to climate variability are selective, with thermal and chemical parameters exhibiting higher sensitivity than dynamic parameters (Insan et al., 2025).
Results of the Kruskal–Wallis test and their interpretation.
| Parameter | p-value | H-value | Interpretation |
|---|---|---|---|
| SST | 0.018 | 8.00 | Significantly different |
| Salinity | 0.024 | 7.45 | Significantly different |
| Nitrate | 0.018 | 8.03 | Significantly different |
| Current velocity | 0.524 | 1.29 | Not significant |
SST, sea surface temperature.
Overall, these results demonstrate that climate variability not only modulates the physical conditions of the ocean surface but also indirectly controls nutrient availability through coastal oceanographic processes. This pattern has important ecological implications, particularly for determining marine productivity and habitat suitability for fisheries and aquaculture activities in the southern coastal waters of Java.
The HSI sensitivity maps illustrate how spatial variations in HSI values occur when one environmental parameter is removed from the model under different climate conditions influenced by the IOD. In general, brighter colors on the maps indicate larger differences in HSI values, meaning that the removed parameter has a stronger influence on habitat suitability. Under Negative IOD conditions, noticeable changes in HSI occur when the nitrate and bathymetry (water depth) parameters are removed, particularly in offshore waters farther from the coastline. In contrast, removing current velocity results in relatively small and spatially homogeneous changes, indicating that current velocity has a relatively minor influence on habitat suitability in the waters of Cilacap and Kebumen during this phase (Figure 10).
HSI Sensitivity Map.
Under Neutral IOD conditions, the sensitivity pattern appears more evenly distributed across the study area. Significant changes in HSI occur when salinity or sea surface salinity (SSS) and bathymetry are excluded from the model, indicating that these parameters play an important role in determining habitat suitability when oceanographic conditions are relatively stable. Meanwhile, during Positive IOD conditions, larger HSI differences appear when salinity and current velocity are removed, particularly in the southern offshore waters. This suggests that water mass dynamics, such as variations in salinity and current transport, become more influential factors during the positive IOD phase.
The spatial patterns observed in the sensitivity maps are consistent with the parameter contribution analysis. The contribution values suggest that bathymetry (water depth) provides a relatively higher contribution compared to other parameters across the examined climate conditions, with contribution values ranging from 0.298 to 0.301, suggesting that water depth may represent an important baseline environmental factor influencing habitat suitability within the context of this HSI model. Salinity and SST also show relatively high contributions, ranging from 0.205 to 0.311 and 0.222–0.26, particularly under Neutral and Positive IOD conditions. In contrast, nitrate and current velocity exhibit relatively smaller contributions, ranging from 0.206 to 0.25 and 0.145–0.278, although they still play a role in shaping habitat suitability variability. Overall, the combination of sensitivity map analysis and parameter contribution values suggests that bathymetry (water depth) may play a relatively important role in shaping habitat suitability patterns, while other oceanographic parameters function as dynamic variables whose influence can vary depending on the prevailing IOD climate conditions (Figure 11). However, these results should be interpreted cautiously because the HSI model applied in this study is relatively simple and primarily exploratory. The contribution values represent relative influences within the model structure rather than definitive ecological causation.
Environmental Parameter Sensitivity of HSI.
The upwelling phenomenon in the southern Java waters tends to reduce the HSI for the cultivation of E. cottonii, particularly during July–September when upwelling intensity increases. Under normal conditions in 2018 and during the positive phase of the IOD in 2019, the strengthening of upwelling was associated with a decline in HSI values compared with June. In contrast, during the negative IOD phase in 2016, when upwelling signals were relatively weak, HSI values remained comparatively stable throughout the July–September period, with dominant values around 0.4. However, this interpretation should be considered with caution because the analysis of negative IOD conditions in this study is based on only a single year of observation, which is insufficient for drawing long-term climatic generalizations. The results of the Kruskal–Wallis test indicate that SST, salinity, and nitrate concentration exhibit significant differences among climate conditions, whereas surface current velocity does not show a significant difference. The combination of sensitivity map analysis and parameter contribution values indicates suggests that bathymetry (water depth) may play a relatively important role in shaping habitat suitability patterns, while other oceanographic parameters (SST, salinity, nitrate, and surface current velocity) function as dynamic variables whose influence can vary depending on the prevailing IOD climate conditions. Nevertheless, the relationship between IOD phases and habitat suitability identified in this study should be interpreted as a preliminary and exploratory indication rather than a generalized climatic pattern. Future studies incorporating longer time series and a greater number of IOD events are necessary to statistically evaluate the consistency of this pattern and to improve the robustness of climate–habitat relationships at broader climatic scales.