Satellite observations of SO₂ have previously been used as a proxy for volcanic eruptions and volcanic ash transport (Thomas and Prata, 2011; Sears et al., 2013). The temporal variability of volcanic plumes produced over the western SPO (10°N–35°S, 125°E–150°W) during the study period of October 2004 to February 2022 was examined using daily level-3 best pixel SO₂ total column products with a 0.25° regular grid, obtained from NASA’s Aura satellite Ozone Monitoring Instrument (OMI), which is available since October 2004 (https://disc.gsfc.nasa.gov/datasets/OMSO2e_003/summary; Li et al., 2017; Li et al., 2020a). The background SO₂ loads (<0.1 Dobson Units (DU); 1 DU = 2.69 × 10¹⁶ molecules cm⁻²) were removed (Li et al., 2020a). To characterize the spatial distribution of volcanic plumes, we also used the level-2 PCA SO₂ total column products (NMSO2-PCA-L2) with a spatial resolution of 50 km, taken from Ozone Mapping and Profiler Suite (OMPS) onboard NASA/NOAA Suomi National Polar-orbiting Partnership (SNPP) satellite, which is available since 2012 (https://disc.gsfc.nasa.gov/datasets/OMPS_NPP_NMSO2_PCA_L2_2/summary). To ensure data quality and maintain pixel numbers, we accepted OMPS SO₂ values with a solar zenith angle less than 75°, cloud cover less than 80%, and SO₂ values higher than 0.1 DU (Yang et al., 2013; Li et al., 2020a). We analyzed the spatial distribution of OMI level 2-PCA SO₂ total column products with a spatial resolution of 13 × 24 km² for the pre-OMPS period (i.e., prior to 2012). These data were re-sampled to 0.25° by 0.25° using MATLAB cubic interpolation methods. The vertical and horizontal transport of ash plumes was detected using level-2 aerosol subtypes data (version 4.x) and level-1 532 nm total attenuated backscatter coefficient measurements (version 4.x) by the Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) onboard the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) satellite (https://www-calipso.larc.nasa.gov; Winker et al., 2009; Vernier et al., 2016; Kim et al., 2018). To track the vertical profile of ash particles from volcanic eruption, we further used daily level-2 aerosol extinction coefficient at 675 nm (version 2) measured by OMPS Limb Profiler sensor on the Suomi-NPP satellite (https://disc.gsfc.nasa.gov/datasets/OMPS_NPP_LP_L2_AER_DAILY_2/summary; Kramarova et al., 2018; Loughman et al., 2018). The aerosol extinction coefficient data have a vertical resolution of approximately 1.8 km. To monitor atmospheric features on volcanic eruptions, the NASA WorldView images of true color corrected reflectance, derived from both SNPP-Visible Infrared Imaging Radiometer Suite (VIIRS) and Aqua-Moderate Resolution Imaging Spectroradiometer (MODIS) were also used (https://worldview.earthdata.nasa.gov/).

The change in phytoplankton biomass over the western SPO was evaluated using Chl-a products, derived using the OC4 blue-green band ratio algorithm, which is applicable for Case 1 waters (i.e., typically oligotrophic and open oceans) (O'Reilly et al., 2000). The Chl-a estimates are 0.25° gridded daily merged products generated by the weighted average method from the GlobColour dataset (https://hermes.acri.fr). To exclude the local impact, we removed the Chl-a values over shallow waters (bathymetry <200 m; Figure S1A ). Furthermore, we estimated ocean PP, using the vertically generalized production model (VGPM) developed by Behrenfeld and Falkowski (1997). The input data required for the PP estimation were obtained from GlobColour datasets for the Chl-a, photosynthetic available radiation, and euphotic depth, and the operational sea surface temperature and sea ice analysis (OSTIA) system dataset for sea surface temperature (SST) (https://marine.copernicus.eu; Good et al., 2020), respectively. The SST data with the spatial resolution of 0.05° were re-sampled to 0.25° by 0.25° grid using MATLAB cubic interpolation methods.

Emitted volcanic ashes are dispersed in the atmosphere and transported by wind to the distance of ten to thousands of kilometers away from their source. The NOAA’s hybrid single-particle Lagrangian integrated trajectory (HYSPLIT) model was applied to identify the atmospheric transport and dispersion of volcanic ash (http://www.ready.noaa.gov) (Stein et al., 2015; Rolph et al., 2017). Global data assimilation system (GDAS) metrological data with a horizontal resolution of 1° were used as model input to calculate the 96-hrs forward trajectories for heights in the surface layer (0.8 km) and stratospheric layer (15 and 21 km) from Ulawun volcano, Papua New Guinea on 7 October, 2006 (4.27°S, 152.20°E), Ambae volcano, Vanuatu on 27 July, 2018 (15.40°S, 167.84°E), Ulawun volcano, Papua New Guinea on 26 June, 2019 (5.05°S, 151.33°E) and Hunga volcano, Tonga on 15 January, 2022 (20.54°S, 175.38°W), respectively. The volcanic ash transport was also examined by using the daily zonal wind (u) and meridional wind (v) at 925 hPa from the National Centers for Environmental Prediction (NCEP)-National Center for Atmospheric Research (NCAR) Reanalysis 1 (https://psl.noaa.gov/data/gridded/data.ncep.reanalysis.pressure.html).

To confirm the influence of the vertical mixing of the water column over the western SPO, we analyzed the weekly estimates of mixed layer depth (MLD) with a 0.25° × 0.25° regular grid derived from the Multi Observation Global Ocean ARMOR3D L4 near real time weekly products and multi-year reprocessed products, for the period of 2020 to 2022 and 2004 to 2021, respectively, distributed by the Copernicus Marine Environment Monitoring Service (CMEMS) (https://resources.marine.copernicus.eu). MLD data were estimated as the minimum value of a density threshold equivalent to a 0.2°C variation of the temperature conditions from the temperature at 10 m depth and the MLD temperature criteria with 0.2°C threshold (de Boyer Montégut et al., 2004). The weekly MLD data were re-sampled to daily data using MATLAB resample function. We also used the daily surface precipitation rate dataset derived from NCEP-NCAR reanalysis 1 to understand the nutrient supply from river flux (Shiozaki et al., 2014). To investigate the nutrient flux from the sediment arounds the islands (Dutheil et al., 2018), we analyzed surface zonal current using Ocean Surface Current Analysis Real-time (OSCAR) dataset with a 1/3 degree grid with a 5 day resolution, which was calculated from satellite datasets using a simplified physical model of an upper ocean turbulent mixed layer (https://podaac.jpl.nasa.gov/dataset/OSCAR_L4_OC_third-deg_YEARLY).

To investigate the role that varying chemical compositions of ash may have had on phytoplankton response, the geochemical composition of volcanic ash from the studied volcanic eruptions was investigated. For the 2018 Ambae eruption, major and trace element geochemistry for the bulk ash and individual ash glasses was taken from Moussallam et al. (2019). Comparable data for the 2009 and 2014–2015 eruptions of Hunga were taken from Brenna et al. (2022), supplemented by analysis of 2022 products collected on land from Tongatapu, Tonga. These samples were collected in the aftermath of the eruption, roughly 75 km south of Hunga volcano. Bulk X-Ray Fluorescence (XRF) analysis of the whole ash sample from Nakualolofa, Tongatapu and XRF analysis of pumice separated from ash collected from Fu’amotu airport was completed by SpectraChem Analytical Ltd. Christchurch using Li-Borate beads for major elements and pressed-powder pellets for trace elements. Individual glass shards from 2022 material deposited on Tongatapu were analyzed using a JEOL Field Emission Electron Probe Microanalyser System 8530F (Hyperprobe) at the University of Auckland. A defocussed beam of 10–20 μm diameter was used with an accelerating voltage of 15 kV, with Na analyzed first and probe conditions monitored using secondary international glass standards.

The spatial distribution of the combined climatology of the OMI-derived mean SO₂ in the western SPO averaged over last 19 years (2004–2022) is shown in Figure S1B . Comparatively higher SO₂ total column values (>0.4 DU) were found in and around the land than that over ocean. The highest values of SO₂ are found in the regions near the volcanoes shown in Figure S1 , indicating the SO₂ variation for the study period in the western SPO is mainly related to the state of volcanic activity. One exception is in Southeast Australia where SO₂ emissions are entirely anthropogenic from 19 coal-fired power stations operating in Australia (Hendryx et al., 2020).

To describe the explosive volcanisms that occurred in the western SPO over last 19 years, we first investigated the volcanic explosivity index (VEI; https://volcano.si.edu/), which distinguishes the volcanic eruptions in the range from 0 (non-explosive eruptions) to 8 (mega-colossal explosive eruptions) on the logarithm of 10, based on volume of magma erupted during an eruption and plume height (Newhall and Self, 1982). During this period, twelve explosive eruptions of VEI 3 or greater have been observed ( Table S1 ). The Hunga eruption from an andesitic submarine caldera, on 15 of January 2022, a Paroxysmal eruption, was the biggest explosion of the 21st century from with a VEI estimated at 5, which is greater by two orders of magnitude than any others in the dataset. However, there were four further events of VEI 4 from the volcanoes in Papua New Guinea, which are Ulawun, Manam, and Rabaul volcanoes. The Ambae event was the most energetic episode of a longer period of sub-Plinian/Plinian volcanism from a basaltic vent at ~1400 m elevation, with VEI of 3 (Moussallam et al., 2019).

Figure 1A shows the daily variations of OMI-derived SO₂ total column during the study period of October 2004 to February 2022 averaged over the western SPO (10°N–35°S, 125°E–150°W). The daily mean SO₂ total column averaged over the western SPO displays generally low values (0.23 ± 0.05 DU), but episodic SO₂ peaks occurred for twelve explosive eruptions of VEI 3 or greater (>0.9 DU) ( Figure 1A and Table S1 ). To identify which explosive eruptions delivered the emitted substances with high intensity over the broad area of the western SPO, we extracted extremely high values fall that above 99.9^th percentile of daily mean SO₂ total column (>0.931 DU), suggestive of high volcanic ash deposition to the western SPO. There were four distinct peaks in October 2006 (maximum value: 1.63 DU), July 2018 (maximum value: 1.56 DU), June 2019 (maximum value: 0.93 DU), and January 2022 (maximum value: 1.95 DU) ( Figure 1B ), corresponding with timing of volcanic eruptions in the Rabaul, Papua New Guinea (October 2006), Ambae, Vanuatu (July 2018), Ulawun, Papua New Guinea (June 2019), and Hunga, Tonga (January 2022). The number of pixels (0.25° grid) with SO₂ values above 99.9^th percentile (0.931 DU) over the western SPO area also showed the strong peaks (>1000 pixels) in the corresponding periods ( Figure S2 ), i.e., covering an area of >765,000 km². The true color images confirm the emissions of volcanic ash from Rabaul (7 October, 2006), Ambae (especially on 20 and 27 July, 2018), Ulawun (25 June, 2019), and Hunga eruptions (14 January, 2022) ( Figure S3 ) (https://worldview.earthdata.nasa.gov). Combined, these data suggest the eruptions of Rabaul, Papua New Guinea (October 2006), Ambae, Vanuatu (July 2018), Ulawun, Papua New Guinea (June 2019), and Hunga, Tonga (January 2022) were the most extreme volcanic plumes produced in the western SPO over last 19 years (2004–2022).

During a volcanic eruption, buoyant plumes of ash and gas are produced (fragmented volcanic glass and minerals, rich in macro- and micro-nutrients; Duggen et al., 2007; Longman et al., 2022). These are rapidly transported horizontally and vertically by the wind, until deposition. Ash deposition in the western SPO may supply limiting nutrients (particularly iron for N₂ fixation; Bonnet et al., 2017) for phytoplankton growth (Heffter and Stunder, 1993; Langmann et al., 2010; Olgun et al., 2013a). To characterize the spatial pattern of the deposition of volcanic materials into the western SPO, in this section, we describe the horizontal and vertical transport and dispersion of ash plumes from Rabaul volcano in October 2006, Ambae volcano in July–August 2018, Ulawun volcano in June–July 2019, and Hunga volcano in January 2022.

Figure 2 shows the horizontal evolution of satellite-derived mean SO₂ total column over the western SPO averaged for the periods of the four volcanic eruptions, 2006 Rabaul eruption (7–13 October), 2018 Ambae eruption (20–27 July and 28 July–4 August), 2019 Ulawun eruption (26 June–4 July), and 2022 Hunga eruption (15–20 January). On 7 October 2006, very high SO₂ columns (>14 DU) were observed around Rabaul volcano ( Figure S4A ), and SO₂ plumes from this source were horizontally distributed over the western SPO, showing two pathways of westward transport and southeastern transport, respectively ( Figures 2A and S4A ). Highest plumes were displayed over near coastal waters with a westward relatively short distance (~700 km) around the northeastern area of Papua New Guinea after leaving the source, during southeasterly surface wind conditions. Following the continuous SO₂ emission from Ambae volcano (>14 DU) during the period of 20–27 July, 2018, with relatively high SO₂ levels around source region (Ambae) (shown in inset in Figure 2B ), the SO₂ plumes averaged during the period of 28 July to 4 August 2018 showed eastward long-distance (>2000 km) transport, showing highest plumes around the northern part of Tongatapu island (29 July 2018), but the direction of surface winds in 2018 case showed different pathways of long-distance transport of SO₂, suggesting that this transport was associated with air masses at the higher altitude ( Figures 2B and S4B ). During 2019 Ulawun eruption event, the SO₂ plumes were relatively less dispersed over the western SPO, compared to those from other volcanoes studied here ( Figures 2C and S4C ). This event showed high plume concentrations only in the eastern area of Papua New Guinea, a short-distance (~500 km) southward from the source region as described in a recent study (McKee et al., 2021). Conversely, the SO₂ distribution averaged during the period of 15–20 January illustrates that the SO₂ emitted from the Hunga volcano travelled a long distance (>6000 km) westward, reaching the northern part of Australia, passing over the waters around Vanuatu and Fiji ( Figures 2D and S4D ). The highest SO₂ values (>14 DU) appeared in around source region (Hunga) and the water surrounded by Vanuatu, New Caledonia, and Fiji (i.e., an area of ~8.2 × 10⁵ km², ~2000 km from source; blue box in Figure 2D ) on 16 January 2022. The transport of volcanic plumes in the 2022 case was consistent with the easterly wind direction at the near surface.

Based on NOAA’ HYSPLIT 96-hrs forward trajectories model, we investigated the vertical transport pathways of tropospheric (0.8 km) and stratospheric air masses (15 and 21 km), originating from the Rabaul, Ambae, Ulawun, and Hunga volcanoes ( Figure 3A ). During the Rabaul volcano event (7 October, 2006), air mass trajectories simulated by HYSPLIT showed that the air mass movement at the altitude of 15 km captured the horizontal distribution of SO₂ plumes, suggesting the stratospheric injection and mass loading of volcanic materials expelled from Rabaul volcano ( Figures 2A , 3A ). The movement of air masses at the altitudes of 0.8 and 21 km from Ambae volcano (27 July, 2018), was consistent with that of near surface wind ( Figure 2B ). Unlike the direction of other air masses, the air mass at an altitude of 15 km showed a pathway consistent with the spatial distribution of SO₂ total column, supporting the suggestion that volcanic ash from the Ambae eruption in July 2018 was transported at an altitude of 15 km (i.e., stratosphere aerosol level). This result is also consistent with the vertical distribution of the Ambae volcanic plume observed with a core brightness temperature (Kloss et al., 2020). The airmass at the altitude of 21 km following the Ulawun eruption (26 June, 2019) showed a pathway eastward, which is consistent with the spatial distribution of SO₂ with low levels ( Figure 3A ). On the other hand, airmass forward trajectory, starting from altitude of 15 km, depicted high SO₂ plume distribution around the eastern area of Papua New Guinea ( Figures 2C , 3A ). The pathways of air masses at the altitudes of 0.8, 15, and 21 km, starting at the Hunga on 15 January 2022, appeared to be in the westward direction, showing consistent pathways to the region of Vanuatu, New Caledonia, and Fiji ( Figures 2D , 3A ). Dispersion of air mass at an altitude of 21 km was nearly the same as that of the SO₂ distribution, indicating that volcanic ash is suspected to be more efficiently transported over a long distance to the northern part of Australia along an altitude of 21 km, higher than the Ambae volcano. Recent studies also showed that much of the plume reached an altitude higher than 30 km during the Hunga eruption, and it developed a massive umbrella region of ~400 km diameter (Carr et al., 2022; Smart, 2022), making it much more effective in dispersing ash over large areas.

To identify whether long-range transportation of air mass in the stratosphere (usually higher than approximately 12 km above the surface) is associated with the dispersion of volcanic aerosol plumes over the western SPO, we further analyzed the CALIPSO L2 aerosol products after the 2006, 2018, 2019 and 2022 eruptions, which classify aerosols into ten subtypes for the stratospheric layer; clean marine, dust, polluted continental/smoke, clean continental, polluted dust, elevated smoke, dusty marine, polar stratospheric clouds aerosol, volcanic ash, sulphate/other ( Figure 3B ) (Kim et al., 2018). The stratospheric aerosols were well spread at the altitudes of around ~15 to 25 km in 2006, 2018, and 2019 eruptions. For the 2022 eruption, relatively higher heights, at the altitudes of 18 to 30 km, respectively, along the CALIPSO measurement orbit tracks, as shown by the HYSPLIT atmospheric long-distance trajectories (15 km in 2006, 2018, and 2019 year and 21 km in 2022 year) and SO₂ distribution ( Figures 2 , 3B ). The presence of volcanic ash plumes was also characterized by CALIPSO L1B total attenuated backscattering coefficient at 532 nm (~0.002 km^-1 sr^-1; yellow in the figure) ( Figure S5 ). Most importantly, the CALIPSO measurements show that the composition of the 2006, 2018, and 2019 plumes was primarily sulphate aerosols (grey in the figures), while the 2022 Hunga eruption was a mixture of volcanic ash and sulphate aerosols (red and grey in the figure), indicating the persistence of ash-laden plumes over the western SPO after the eruption ( Figure 3B ). Furthermore, we investigated the vertical distribution of OMPS L2 aerosol extinction to understand the quantitative retrievals of volcanic ash transported and deposited to the western SPO following 2018, 2019, and 2022 eruptions ( Figure 3C ). There was no available data in 2006. The vertical distribution of aerosol extinction along the OMPS orbit track on 29 July 2018, confirmed a slight increase (~0.02 km^-1) at altitudes of ~15 to 20 km in the northern area of Tongatapu island where the highest levels in SO₂ were shown, indicating relatively potential low deposition ( Figures 2B , 3C ). After the 2019 volcanic eruption, aerosol extinction coefficients on 27 June 2019 also showed a slight increase (>0.02 km^-1) at altitudes of ~15 to 20 km around Papua New Guinea. However, unlike these events, the vertical distribution of aerosol extinction coefficients on 16 January 2022 showed relatively higher values with levels of >~0.05 km^-1 at all observed altitudes (<30 km) in the water surrounded by Vanuatu, New Caledonia, and Fiji ( Figure 3C ), where ash-laden plumes and high levels in SO₂ were mainly observed, but relatively low aerosol extinction coefficients in the western part of Vanuatu on 17 January 2022 ( Figures 2D , 3A, C ). Mishra et al. (2022) showed that on 16 January the maximum SO₂ value in the region occurred when the spread speed of atmospheric westward plume was lowest as <1000 km day^-1, but >2500 km day^-1 after 17 January. Ash deposition maxima can result from the trapping of ash related to relatively low wind speed, as shown in a numerical modeling study (Poulidis et al., 2018). Therefore, these results suggest 2022 Hunga eruption was accompanied by high deposition of ash-rich aerosols over the water surrounded by Vanuatu, New Caledonia, and Fiji on 16 January 2022, probably by relatively weak winds.

The composition of the deposited volcanic ash has the potential to control the phytoplankton response, with differing levels of nutrient release from different ash compositions already noted (Jones and Gislason, 2008; Longman et al., 2022). Our ash classification showed that in October 2006, the large sub-Plinian eruption from Rabaul volcano, with a whole-rock composition of Trachydacite ( Figure S6 ) (Bouvet de Maisonneuve et al., 2015; Bernard and Bouvet de Maisonneuve, 2020). The 2018 Amabe volcano, which is a basaltic shield volcano, produced basaltic andesite to trachy-basalt ash (Moussallam et al., 2019). The 2022 Hunga eruption produced andesitic to basaltic andesite ash, similar to earlier eruptions of this volcano (Brenna et al., 2022). Figure S7 shows that both Ambae and Vanuatu ash contain high iron contents in their total glass and bulk analyses (often >10 wt%), which is a common nutrient deficiency in the ocean, but Rabaul ash contains relatively low iron contents (<10 wt%) ( Table S2 and S3 ) (Moore et al., 2013). This may explain the lack of phytoplankton response to Rabaul, but our data do not contain information on bioavailability of the Fe contained within the ash. As such, further research is necessary to confirm the hypothesis that changing ash composition leads to different phytoplankton response.

To understand the response of phytoplankton subsequent to distal ash transport from four extreme eruptions (2006 Rabaul, 2018 Ambae, 2019 Ulawun, and 2022 Hunga eruptions) in the western SPO, we carried out an inspection of area-averaged time series of GlobColour-merged Chl-a concentrations in the distal region ( Figures 4A–D ), i.e., regions having highest atmospheric SO₂ loading (>14 DU) extending toward the western SPO (blue box; 2006: 0–4°S, 145–150°E, 2018: 12–18°S, 178°E–169°W; 2019: 6.5–10°S, 147–157°E; 2022: 15–23°S, 165–174°E) ( Figure 2 ). A time series of Chl-a estimates in the distal region, which is affected regions by long-range transport of volcanic plumes (>500 km), showed profoundly different response to 2022 Hunga events when compared to all other volcanic events (2006, 2018, 2019 events; Figures 4A–D ). Following long-range transport of the 2006, 2018, 2019 volcanic emissions ( Figures 2A–C ), Chl-a concentrations did not show any noticeable increases ( Figures 4A–C ), suggesting limited nutrient supply (particularly iron) probably due to low volumes of ash deposition. However, in regions affected by likely high deposition of long-distance transported Hunga ash on 16 January 2022 ( Figures 2D , 3C ), Chl-a estimates were two times higher than that of climatological mean and one standard deviation for approximately 10 days after a lag period of 6 days ( Figure 4D ). The difference between Chl-a values from the 2022 Hunga event with the climatological values were statistically significant (p <0.05; t-test analysis). The response of Chl-a values to 2022 Hunga eruptions (0.15 mg m^-3) showed a 67% increase in average compared to climatological values (0.09 ± 0.03 mg m^-3) (red shading in figure). Oceanic PP estimates subsequent to 2022 Hunga eruptions (383 mg C m^-2 d^-1) also showed a 43% increase compared to climatological value (268 ± 74 mg C m^-2 d^-1). This conclusion is echoed by the findings of a recent study, which showed increasing Chl-a in close proximity to Hunga volcano, but with a higher magnitude in Chl-a values (>10-fold than before eruption) immediately after the eruption (Barone et al., 2022). However, as the presence of ash particles can directly influence the satellite ocean color signatures, in turn biasing the appraised Chl-a estimates, monitoring of the phytoplankton response to the ash deposition using satellite-derived optical datasets should be interpreted with caution (Johnson et al., 2011; Browning et al., 2015). The ash particles have residence times of up to 2 days in the euphotic zone depending on the particle sizes (Duggen et al., 2007). Whiteside et al. (2023) revealed that there were extremely high amounts of suspended ash particles in the surface waters around the Hunga volcano 2 days after 2022 eruption, which probably led to ash-contaminated Chl-a estimates but 9 days after the eruption the ash-related values recovered to normal, suggesting the removal of ash signal in the surface with time. In addition, previous studies have shown that the phytoplankton began to respond five to six days after deposition of volcanic ash, as shown in mesoscale iron-fertilization experiments (Coale et al., 2004; Duggen et al., 2007; Langmann et al., 2010; Lin et al., 2011; Yoon et al., 2018). Therefore, in the context of this study, the increase in Chl-a at a time period ∼6 days after 2022 Hunga eruption, potentially after sinking of ash particles below surface waters, would exclude the misinterpretation of Chl-a estimates by ash particles itself in the surface water (Langmann et al., 2010).

Furthermore, we analyzed the spatial distribution of the difference of 10-days mean Chl-a concentrations 6 days after volcanic eruption (2006: Chl-a_{2006 (13–22 October)}; 2018: Chl-a_{2018 (2-11 August)}; 2019: Chl-a_{2019 (2–11 July)}; 2022: Chl-a_{2022 (21–30 January)}) and 10-days mean climatological Chl-a concentrations (2006: Chl-a_{2006clim (13–22 October)}; 2018: Chl-a_{2018clim (2-11 August)}; 2019: Chl-a_{2019clim (2–11 July)}; 2022: Chl-a_{2022clim (21–30 January)}) averaged for corresponding periods for 19 years (2004–2022) except corresponding year, respectively (e.g., ΔChl-a _{2006 (13–22 October)} = Chl-a_{2006 (13–22 October)} – Chl-a_{clim (13–22 October) except 2006}) ( Figures 4E–H ). As shown in the temporal analysis, the long-distance transport of 2006, 2018, 2019 volcanic plumes over western SPO was not followed by any apparent response in Chl-a over regions that have high levels in SO₂ ( Figures 2A–C , 3C , 4E–G ). However, in 2022, considerable positive anomalies of Chl-a values ( Figure 4H ) were also distinctively visible over the water surrounded by Vanuatu, New Caledonia, and Fiji which showed potential high deposition following the long-range transport of ash released from Hunga ( Figures 2D , 3C ). Furthermore, an anomalously high phytoplankton stock was detected in waters in close proximity to Hunga volcano at a time period 6 days after the 2022 eruption ( Figure 4H ), as shown in Barone et al. (2022), indicating a continued response in Chl-a even under a substantial sinking of ash particles (Duggen et al., 2007; Whiteside et al., 2023). These results suggest that volcanic ash deposition is considered the likely explanation for the enhancements in phytoplankton response in the western SPO.

In the western SPO, PP has been linked to the occurrence of N₂ fixing organisms, and high N₂ fixation rates have been recently assumed to be due to the alleviation of iron limitation, driven by multiple potential iron inputs, such as shallow hydrothermal plumes, island sediment, and river runoff (Shiozaki et al., 2014; Caffin et al., 2018; Dutheil et al., 2018; Guieu et al., 2018; Tilliette et al., 2022). Therefore, to investigate the hypothesis that significant ash deposition on the western SPO led to the Chl-a increase following the 2022 Hunga eruption, we investigated the possibility of Chl-a response to nutrient supply from other sources.

First, the shallow hydrothermal sources along Tonga-Kermadec arc have the potential to fertilise the western SPO with iron (Guieu et al., 2018; Tilliette et al., 2022). However, as Tilliette et al. (2022) say, shallow hydrothermal plumes release high levels of iron only around the shallow hydrothermal sources on the Tonga-Kermadec arc, and to lesser extent, into the water surrounded by Vanuatu, New Caledonia, and Fiji. This indicates the Chl-a increase following the 2022 Hunga eruption was not driven by a shallow hydrothermal source of iron.

Dutheil et al. (2018) showed that the spatial distribution of Chl-a estimates is tightly controlled by iron release from the sediment flux, which are related to the zonal (eastward) advection of iron downstream of the islands, a second potential iron source. Figure S8A shows the spatial distribution of the 5-day mean surface zonal current anomalies for January 2022 against 2004–2021 climatology. Positive (negative) values of surface zonal currents imply eastward (westward) anomalies. However, the surface zonal current dataset showed the insignificant changes in this region between pre- and post-eruption of Hung volcano in January 2022, when compared to other years.

Many islands could also deliver terrigenous nutrients through runoff to the western SPO. A large amount of runoff from land is expected to be triggered by high precipitation, which causes the water column stratification and land drainage required for N₂ fixation (Shiozaki et al., 2014). Figure S8B shows the daily variation of precipitation rates averaged over the western SPO for 2022 vs. climatological mean with one standard deviation of the mean (2004 to 2021). However, there was no significant increase in precipitation over the waters of Vanuatu, New Caledonia, and Fiji during 2022 volcanic episodes compared to climatological values, indicating a likely minor impact of river runoff on nutrient supply.

A final alternative explanation for the changes in phytoplankton biomass in the western SPO may be the changes in MLD, related to the supply of nutrients from below the thermocline to the euphotic zone (Vaughan et al., 2007; Boyce et al., 2010; Mantas et al., 2011; Jutzeler et al., 2014; Jutzeler et al., 2020). Indeed, Terry et al. (2022) reported the tsunamis produced by the 2022 Hunga eruption, were up to 15 m on the islands near to the source volcano. However, there was no significant change in the area-averaged MLDs in 2022 year with time, relative to climatological area-averaged MLDs ( Figure S8C ). Therefore, in summary, the lack of alternative plausible hypotheses supports our conclusion that the increase in Chl-a is linked to the 2022 Hunga ash fall over the western SPO.