Field data from multiple sites spanning the WAP were utilized in this study ( Figure 1A ). The majority of our combined field dataset was from the Palmer Long Term Ecological Research (LTER) grid stations ( Figure 1B ). Discrete δ¹⁸O samples were collected from oceanographic cruises, primarily in austral summer (December to February), conducted between 2011 and 2019 over the WAP shelf. Samples were drawn from Niskin bottles closed at discrete depths during water column profiling with a Sea-Bird conductivity–temperature–depth (CTD) instrument (Ducklow et al., 2013; Meredith et al., 2017a).

In addition to water sample collection over the WAP shelf, the Palmer LTER program also samples quasi-weekly during austral summer from two nearshore sites near the Palmer Research Station on Anvers Island ( Figure 1B ). Station B, located at 64.7795°S/64.0725°W, is <1 km from shore and in waters of 75 m depth. Station E, located at 64.815°S/64.0405°W, is farther offshore and more directly exposed to the continental shelf and open ocean, with a water depth of 200 m. The start of each sampling season is subject to boat access due to severe weather conditions and sea ice presence, thus ranging from mid-October to late December. During austral winter, when boat operation ceases, the Palmer LTER program also samples seawater through an inlet at 5.8 m depth, located at 64.7738°S/64.0545°W. These nearshore stations are situated near the shallow marine-terminating Marr Ice Piedmont and other land-based ice masses (Meredith et al., 2021) ( Figure 1B ).

The combined field dataset also contains δ¹⁸O seawater samples from additional nearshore stations at Potter Cove and Maxwell Bay, located on King George Island ( Figure 1B ). Water samples were collected using 4.7 L Niskin bottles during water column profiling with a Sea-Bird SBE 19 CTD instrument (Meredith et al., 2018). Some field data from the Rothera Oceanographic and Biological Time Series (RaTS) were also included in the combined field dataset, albeit sparse due to a lack of matching MODIS-to-field data given its proximity to the shore of Adelaide Island (67.5700°S/68.2250°W) ( Figure 1B ). Water sampling here was also conducted in conjunction with CTD casts (Meredith et al., 2010). This expansive selection of study sites aligns with prior studies, such as the WAP research overviews provided by Ducklow et al. (2007) and Hendry et al. (2018), which often include islands and areas near the northern tip of the Peninsula.

From each of the sites mentioned, δ¹⁸O seawater samples at the surface were selected for this study. Surface is here defined as the shallowest sample depth in the water column within the top 5m, and δ¹⁸O is here measured in units of ‰. SSGM is derived based on mass balance calculations conducted in the aforementioned studies, and it is measured in units of freshwater % (Meredith et al., 2010; Meredith et al., 2017a; Meredith et al., 2018; Meredith et al., 2021). Each discrete δ¹⁸O sample is paired with its corresponding salinity values from CTD data. The mass balance calculation presumes each sample is composed of a simple mixture of three components – ocean water (ow), sea ice meltwater (sim), and meteoric water (met), with the latter term being the sum of precipitation and glacial meltwater:

This system of equations is solved for F_sim , F_met , and F_ow , which are the respective fractions of the three components in each sample. S_sim , S_met , and S_ow are the salinity values for the end-member source components, while δ¹⁸O_sim , δ¹⁸O_met , δ¹⁸O_ow are the corresponding δ¹⁸O values. Endmember values are taken from the previous studies that generated our combined field dataset ( Table 1 ); these data sources include Palmer LTER grid (Meredith et al., 2017a), Palmer Station on Anvers Island (Meredith et al., 2021), Potter Cove (Meredith et al., 2018), and Rothera Point on Adelaide Island (Meredith et al., 2010) as described above.

MODIS-Aqua Level 3 remote sensing reflectance data, R_rs(λ), for 10 spectral bands (i.e., 412, 443, 469, 488, 531, 547, 555, 645, 667, 678 nm where λ is light wavelength in vacuum) were acquired from the NASA Ocean Biology Distributed Active Archive Center (DOI: 10.5067/AQUA/MODIS/L3M/RRS/2018). Level 3 Ocean Color Standard Mapped Images at 4 km resolution were accessed via Google Earth Engine (GEE). An initial matchup dataset was compiled by iteratively reviewing all field data, and then matching each field data point with a set of corresponding MODIS pixel values that cover these field coordinates from the same date. Furthermore, the level 3 R_rs(λ) data already removes areas with sea ice, thus eliminating the potential of sea ice contamination in the overall dataset. This process produced a dataset where each field measurement of meteoric water fraction was matched with its ten corresponding MODIS R_rs band values. Out of the 86 matchups, 38 were from the Palmer LTER grid, 42 from the Palmer Station site, and 6 from Potter Cove or the Rothera Time Series. The low matchup availability is mostly due to cloud cover. This initial dataset also needs very low or zero glacial meltwater values for standard machine learning models to effectively detect a wide range of target values. Therefore, the initial dataset was also supplemented with synthetic “blank” meteoric water samples. For each date when a field-MODIS data point pair is found, a synthetic blank sample is also obtained at three locations on the outside of the Antarctic Circumpolar Current (ACC) ( Figure 1A ). They were sampled at 47.2382°S/99.1695°W, 54.3690°S/97.1654°W, and 35.5635°S/129.8333°W. For each “blank” sample, a value of 0% glacial meltwater fraction was assigned to the R_rs band values obtained for these pixels.

In this study, the MODIS R_rs(λ) band values are input variables while meteoric water fraction is the target variable. The input variables are the ten ocean color band values of R_rs(412), R_rs(443), R_rs(469), R_rs(488), R_rs(531), R_rs(547), R_rs(555), R_rs(645), R_rs(667), R_rs(678), and band math values for R_rs(555) + R_rs(667), and R_rs(667)/R_rs(488) per previous AOP models for retrieving glacial meltwater fraction (Pan et al., 2019). In order to model the data, we executed a three-step process involving: (1) algorithm selection, (2) hyperparameter optimization of the selected algorithm using train, validation and test subsets, and (3) final model training with the best hyperparameters on the complete dataset. For algorithm selection, we used a method of evaluating thirty-six different regression algorithms, including adaptive boosting, support vector machines, k-nearest neighbors, linear regression, gradient boosting, an individual decision tree and a random forest, among others (https://scikit-learn.org/stable/supervised_learning.html). The selected algorithm had to be reasonably interpretable, so algorithms like advanced neural networks were excluded. The primary goal of this preliminary step is to understand how different models would respond to the SSGM training dataset as well as their general performance without extensive model tuning, such as hyperparameter search.

After identifying the random forest algorithm as the best-performing in the initial evaluation, we conducted a hyperparameter search to fine-tune the following hyperparameters: n-estimator, which is the number of decision trees in a random forest; minimum-samples-leaf, which is the minimum number of samples that should be present in the leaf node after splitting a node; and minimum-samples-split, which is the minimum number of data points in any given node in order to split it.

Hyperparameter optimization is the process of iteratively tuning a learning algorithm’s hyperparameters in order to better fit training data. Overfitting is a common problem when utilizing this technique (Cawley and Talbot, 2010; Feurer and Hutter, 2019). To avoid overfitting, the dataset was split into training (80%), validation (10%) and test subsets (10%) ( Figure 2 ). The training subset was used to train 10000 candidate models. The cross-validation accuracy score for each model was evaluated using the validation dataset. There are a number of search algorithms used in hyperparameter optimization – most commonly grid search, random search and Bayesian algorithms. A Tree-structured Parzen Estimator (TPE), which is a Bayesian technique, enables a more efficient search of the hyperparameter space than grid or random searches (Bergstra et al., 2011). This means new hyperparameters were selected iteratively based on the previous performance of those hyperparameters, as measured against the validation dataset. The best hyperparameters were 265 n-estimators, a minimum-samples-leaf of 1, and a minimum-samples-split of 8. The maximum depth of the random forest is not predetermined, and the decision trees’ nodes are expanded until all leaves contain less than the minimum-samples-split. At the completion of the hyperparameter search, a model was selected and evaluated against the test dataset. This means that the test dataset was used only once during the entire modeling process ( Figure 2 ).

The finalized model is trained using data and settings described earlier in this section – all available R_rs bands and band math as well as the random forest model hyperparameters determined from the TPE-based search process. The finalized model is deployed with geemap (Wu, 2020), a Python package for interactive mapping with GEE, using the native GEE random forest function. The MODIS-Aqua level 3 dataset was compiled in GEE using its NASA/OCEANDATA/MODIS-Aqua/L3SMI image collection (DOI: 10.5067/AQUA/MODIS/L3M/RRS/2018). Further model validation was conducted to assess the model’s robustness and to gain insights on its applicability to spatially and temporally averaged MODIS data ( Figure 3 ). Spatially, each field data point is compared with the MODIS-based prediction that encompasses the field data coordinates ( Figure 3 row 1, column 1), and then the field data point is also compared with spatially averaged MODIS-based predictions – by adding 2.5 km, 5 km, 7.5 km, and 10 km to each side of the pixel to create a new area of interest, and then averaging all predicted values that fall within the new target areas respectively ( Figure 3 rows). Temporally, each field data point is compared with the MODIS-based prediction from the same date ( Figure 3 row 1, column 1), and then the field data point is also compared with predicted values averaged over 3 days, 7 days, monthly, and date range covering the Palmer LTER cruises which varies annually ( Figure 3 columns). These comparisons generate a matrix of 25 configurations of different spatial and temporal correlations between field data and predicted values ( Figure 3 ). For each comparison configuration, we evaluated the Pearson’s correlation coefficient (r), p-value, sample size and a visualization of the correlations as scatter plots ( Figure S1 ). Additional statistical measures were produced to evaluate the SSGM model by focusing on the correlation between the daily MODIS-based values and discrete field data ( Figure 3 row 1, column 1). The loss of correlations between field data and spatiotemporally averaged predictions can provide some indications on SSGM’s spatial scale and residence time.

The final SSGM model was also applied to airborne imaging spectroscopy data from the NASA/JPL Portal Remote Imaging SpectroMeter (PRISM) (Mouroulis et al., 2014), collected as part of the O₂/N₂ Ratio and CO₂ Airborne Southern Ocean (ORCAS) campaign in January-February 2016 (Stephens et al., 2018), to further assess the model’s applicability across platforms. The goal was to demonstrate the spatial resolution advantages of such airborne datasets. PRISM flew on the NSF/NCAR Gulfstream-V aircraft at an approximate altitude of ~12.2 km providing water-leaving radiance and remote sensing reflectance from 350 to 1050 nm, approximately every 3 nm, with 10 m spatial resolution (https://prism.jpl.nasa.gov/dataportal/). Atmospheric correction was conducted using an Optimal Estimation formulation that simultaneously models surface and atmospheric reflectance contributions from statistical priors (Thompson et al., 2018; Thompson et al., 2019). The PRISM ORCAS data were pre-processed in ENVI^® to extract R_rs(λ) bands according to MODIS wavelengths and the images were exported to GEE as GeoTIFF files. Here, we assume that the PRISM bands near the targeted MODIS bands would perform similarly. The final SSGM model was applied to each airborne image to retrieve glacial meltwater fraction.

When the predicted glacial meltwater fraction values and field data are compared during the development phase ( Figure 2 ), we found the validation dataset had an r value of 0.878 and the test dataset had an r value of 0.885. The r values across both datasets are relatively high, and notably, these r values are relatively consistent indicating the hyperparameter search did not cause the model to overfit against the validation data. This evaluation confirms the model’s overall performance, so all development datasets were combined to produce the final model for this study. After the final version of the model is applied to MODIS-Aqua data, we found relatively high r value when field data are compared to predictions based on matching pixel values from the same dates (r = 0.82, p < 0.01). However, the correlation decreases as the MODIS-based predictions are averaged over longer time periods ( Table 2 , columns). The decrease in correlation is even more pronounced when the predicted values are averaged over larger areas ( Table 2 , rows). Figure S1 presents a visualization of these correlations as scatter plots.

To further evaluate the model performance, several statistical measures were calculated to describe the regression function between field data and daily MODIS-derived SSGM values (without spatial or temporal averaging, Figure S1 row 1, column 1). These measures are commonly used for evaluation in biogeochemical modeling and bio-optical studies (Brewin et al., 2015; Seegers et al., 2018) ( Figure 4 ; Table S.1 ). We found that the residual differences between MODIS-derived SSGM and field data have a symmetric probability distribution around a median of 0.04%. The median ratio (MdR) and median bias (MdB, in units of predicted variable) measures are produced to characterize systematic deviations, while the median absolute percentage difference (MdAPD) and root mean square error (RMSE, in units of predicted variable) are used to characterize random deviations between modeled SSGM values and field data. A model that provides the best overall fit to in-situ data should have MdR close to 1 and low values of MdB, MdAPD and RMSE. We found MdR = 0.990, while MdB = −0.040, MdAPD = 6.375%, and RMSE = 0.463. Moreover, errors in remote sensing and bio-optical models are also often reported by statistical indicators calculated based on linear regression between in-situ and predicted values. Here, we calculated the statistical measures based on Model II least-squares fit (Ricker, 1973), because both satellite- and δ¹⁸O-derived values presented in Figure 4 are collected from an uncontrolled environment, subject to error, and are dependent on the actual SSGM fractions in the field. The slope of the line is determined by calculating the geometric mean of the slopes from the regression of Y-on-X and X-on-Y (Sokal and Rohlf, 1995). The best-fit coefficients of the regression function and correlation coefficient reveal the degree to which the predicted values agree with measured values over the entire dynamic range. For instance, deviations from 1 for the slope and deviations from 0 for the intercept (A) of the fitted regression reveal the variation in the bias of SSGM predictions relative to the field data across its full range ( Table S.1 ). These measures are further evaluated when our analysis considers outliers. The outliers were selected based on a Bland-Altmann plot ( Figure S2 ) and they are defined as data points where the predicted SSGM value is beyond ±2 standard deviations from the mean of residuals between predicted and observed values.

Further model diagnostics were conducted to evaluate the model’s performance and to better understand the relationship between the R_rs(λ) values and predicted SSGM fraction values. In Figure 5 , the model’s input variables are ranked by the amount of variance that they explain. The top variable, R_rs(667)/R_rs(488), explains 82.96% of total variance. The following two variables – R_rs(443) and R_rs(412) – explain an additional 7.35% cumulatively. Therefore, the top three input variables explain a total of 90.31% of the variance. The remaining input variables explain 9.69% of the total variance ( Figure 5 ).

Notably, we intentionally introduce R_rs bands that were designed for terrestrial and atmospheric applications – in particular, Band 1 centered around 645 nm, Band 3 centered around 469 nm, and Band 4 centered around 555 nm. The results demonstrate the model’s ability to distinguish valid information (i.e., various R_rs bands designed for ocean observation were selected as top input variables and they can explain 90.31% of the model’s total variance). These results also highlight the model’s robustness to provide valuable insights instead of merely over-fitting to the target variable with all available input data indiscriminately.

The response of each input variable is also examined in more detail. We use individual conditional expectation plots to present the responses of decision trees within the random forest model to changes in each input variable (Molnar, 2020). In each panel of Figure 6 , the blue lines represent how individual decision trees predicted SSGM fraction values change according to the input R_rs(λ) values, while the red line represents the average prediction of these decision trees ( Figure 6 ). Within the top five input variables ( Figure 5 ), predicted SSGM values increase abruptly when R_rs(667)/R_rs(488) is over 0.05 meanwhile a more gradual decrease in SSGM is observed for input variables R_rs(443), R_rs(412) and R_rs(469) when they increase from 0.005 sr^-1 to 0.01 sr^-1. Additionally, a gradual increase in SSGM is observed when R_rs(645) increases from 0.00025 sr^-1 to 0.0005 sr^-1. In contrast, the model exhibits minimal change in predicted SSGM values according to changes in some input variables, such as R_rs(678), R_rs(667), and R_rs(555) ( Figure 6 ).

The final model is applied to MODIS-Aqua data to provide visualizations of SSGM. Figure 7 represents the December, January, and February monthly climatology between 2010 and 2020. The austral summer months coincide with most of the sampling time periods of our field data. SSGM is most abundant along the coast and extends across the WAP continental shelf into the open Southern Ocean. High SSGM is also observed around islands and nearshore locations around outlets of glacial fjords. Notably, some high SSGM is also observed near the Larsen C Ice Shelf on the eastern side of the Peninsula. Overall, there appears to be more glacial meltwater on the southern end of the WAP, particularly in Marguerite Bay and near the Wordie Ice Shelf ( Figures 7A, D, G ). The SSGM monthly climatology is also compared to surface chl-a concentration and sea surface temperature (SST) from the same time periods. High SSGM values appear to coincide with chl-a concentrations, but there is a distinct difference in their spatial patterns over the entire WAP ( Figures 7B, E, H ). Similarly, some nearshore high SSGM is found to coincide with low SST, especially in Marguerite Bay and between Adelaide Island and Lavoisier Island ( Figures 7C, F, I ).

The SSGM model was also applied to PRISM airborne ocean color data to assess the model’s applicability across remote sensing platforms. PRISM scene prm20160125t181722 is presented in Figure 8 , collected NW of Anvers Island. The entire scene has an average of 3.94% SSGM fraction. A few icebergs < 1km in size can be found in the middle of the scene ( Figure 8C ) and, in contrast, the maximum SSGM fraction reaches 6.49% in the vicinity of the icebergs ( Figure 8D ). Overall, the SSGM derived from airborne data shows consistent spatial distribution as the MODIS-based values. In particular, a distinct front can be observed in both datasets immediately south of 64 °S ( Figure 8A ). This consistency in the SSGM spatial distributions derived from different remote sensing platforms illustrates the potential for applying the model to other ocean color datasets beyond MODIS. Moreover, the airborne data demonstrates its significant advantage of having finer spatial resolution; for example, meltwater from the icebergs cannot be observed in the spatially coarse MODIS-based SSGM data, but it can be detected using the airborne data. Notably, a trail of bergy bits (Dowdeswell and Forsberg, 1992) from the icebergs can be clearly seen in both enhanced RGB image and SSGM data from PRISM ( Figures 8C, D ). There are also some discrepancies between the PRISM- and MODIS-derived SSGM values. The differences are likely due to the different processing methods of the remote sensing data from the two different platforms (Material & Methods 2.3, 2.4, and 2.5), The discrepancies are also likely resulted from their temporal difference between the times of data acquisition, which is several hours apart.

The final model is applied to MODIS-based predictions over different spatial and temporal scales. A notable feature is the formation of a “threshold” that is parallel to the x-axes when predicted SSGM values are averaged over extended spatial and temporal scales ( Figure S1 ). This horizontal cluster of data points indicates that it is unsuitable to compare a single field data point with averaged SSGM prediction over extended space and time periods. Because of this, future use of this model should only involve other datasets that are on a similar spatial and temporal scale as the SSGM data product. For an instance, daily 1 km resolution SSGM product should not be directly merged with 8-day 4 km resolution chl-a product; a SSGM composite should be created first before merging the datasets. The rapid decrease in r values across the comparison matrix ( Table 2 ) also represents a potential indication of residence time and spatial extent of SSGM features. While this speculation can only be confirmed with additional studies beyond the scope of this project, it still offers some information on the spatial and temporal requirements for future field campaigns and remote sensing missions to effectively resolve SSGM features.

When the final model’s predictions are compared with the field data, most predicted values appear to slightly over-estimate SSGM, but within an acceptable range ( Figure 4 ). In contrast, outliers appear to significantly underestimate SSGM (except for one sample, Figure S2 ). Statistical measures shown in Table S.1 have been carried out to analyze the dataset, with and without the outliers, and both analyses show relatively consistent results ( Table S.1 ; Figure S4 ). The slope of the regression line between modeled and in-situ values decreased from 0.931 to 0.771 after outliers were removed. This indicates that the model exhibits some remaining patterns in bias over the entire range of data, as seen in Figure S2 and Figure S4 . They may be related to potentially higher measurement uncertainties in some data ranges (e.g., low values of SSGM) or result from uncertainties in spatio-temporal matching of satellite and in-situ measurements. We note, however, that apart from these regression line metrics, all other statistical indicators in Table S.1 of aggregate model bias and random error improved following the removal of these outliers, and overall support the potential utility of this preliminary model developed based on a relatively limited training dataset. Additional observations, particularly in the lower and upper ranges of SSGM, are needed to better assess model uncertainties and improve model performance across the entire SSGM range.

There are inevitably some assumptions and uncertainties associated with the SSGM model. One area of uncertainty is the quantity of the final (combined) training dataset (n = 204, including synthetic blanks). While this dataset is relatively small in comparison to machine learning projects in industry, the range of field data values in the dataset helps train a representative and robust model. Within the dataset, high meteoric water content was captured near Potter Cove. Many of these surface meteoric water samples were comprised of almost entirely glacial meltwater associated with heavy sediment discharge at this site (Meredith et al., 2018). These high meteoric water contents during the growing season are also consistent with the low salinity and high turbidity observed in the nearby Admiralty Bay (Osińska et al., 2023). Low SSGM content was introduced to the dataset via synthetic blanks that were generated outside of the ACC in the Southern Pacific open ocean ( Figure 1 ). This gradient of SSGM captures the full range of meltwater content and their corresponding remote sensing reflectance properties. Moreover, the training dataset in this study is one of the most comprehensive δ¹⁸O data compilation for the WAP region, including four long-term time series at different locations (Material and Methods).

There are also assumptions and uncertainties associated with the field dataset itself. Because the meteoric water content was computed through mass balance calculations (Material and Methods), these values are sensitive to the selection of end member values (i.e., pure glacial and sea ice meltwater δ¹⁸O and salinities). To mitigate these uncertainties, a consistent and rigorous methodology was used during these field studies to select end members that were locally or regionally representative (Meredith et al., 2017b). In addition, the meteoric water variable is comprised of precipitation and glacial meltwater (Material and Methods). While the current state of this method cannot partition these two components directly, δ¹⁸O sampling remains one of the most cost-effective and accessible methods for estimating meltwater content and it can be collected on most oceanographic and community science cruises (Cusick et al., 2020). Moreover, the meteoric water variable represents an upper limit of glacial meltwater estimation, which is still useful for obtaining insights on its surface distribution over space and in time.

Additionally, the “synthetic blanks” in this study (Material and Methods) were selected outside of the ACC, therefore their meteoric water content likely include no glacial meltwater but various amounts of precipitation. Since these blank samples have significantly lower optical signal in comparison to the in-situ data collected at the WAP ( Figure S3 ), by setting the synthetic blanks as 0% SSGM, we have likely trained the model to neglect the effect of precipitation on SSGM predictions. SSGM prediction also appears to be only weakly sensitive to the selection of blank SSGM data ( Figure S3 ). This assumption can be verified to some degree when daily precipitation data is compared with in-situ meteoric water fractions around Anvers Island ( Figure S4 ). Precipitation data retrieved between 1989 and 2019 shows high meteoric water content is associated with low precipitation ( Figure S4A ). Furthermore, the majority of this dataset is associated with no precipitation. Out of 410 data points, 219 are associated with 0 mm of rainfall, 361 are associated with 0 cm of snow precipitation, and 282 are associated with no accumulation by snow stake ( Figure S4A ). In addition, the anomalous values that fall away from the validation line in Figure 4 have been identified as outliers, and these values are also associated with low precipitation ( Figure S4B , S4C ). Hence any predictions that deviate from in-situ data is likely not due to excess precipitation interfering with glacial meltwater signals in δ¹⁸O samples. Due to the complex relationship between precipitation and the δ¹⁸O measurement, such as land-accumulated snow melt discharging into the ocean during summer or any precipitation events that might have occurred shortly before in-situ sampling, data presented in Figure S4 do not definitively decouple glacial meltwater from precipitation. However, these additional results indicate that δ¹⁸O has value as an effective in-situ measurement for estimating glacial meltwater in the WAP, and for training models for predicting SSGM. In the future, additional cost-effective tracers should be explored to enable partitioning of meteoric water content.

The SSGM model has provided important results with broader implications. For instance, the variable importance ranking shows the model is primarily dependent on a few key R_rs wavelengths, indicating these wavelengths are important to capturing the underlying optical characteristics of SSGM ( Figure 5 ), though to a lesser extent, additional bands can also contribute to an increase in overall prediction power of the model, illustrating the merit of utilizing hyperspectral airborne ocean color data for future SSGM model development. Hyperspectral airborne data have several other advantages – such as its high spatial and temporal resolutions which reveal detailed features of SSGM distribution that cannot be resolved with MODIS data ( Figure 8 ), as well as its increased spectral sampling that could provide additional spectral features not seen in multispectral data that can further tease out SSGM from other optical constituents. These results highlight the benefit of supplementing existing multispectral spaceborne datasets with hyperspectral airborne data in order to achieve better detection and monitoring of SSGM. Moreover, the application of this model on PRISM imagery indicates the applicability for other ocean color data products as well; in general, the model can generate SSGM results if the input R_rs data are at or near the MODIS bands ( Figure 5 ) and are processed with consistent techniques used to produce MODIS Level 3 data products. Furthermore, the applicability of this SSGM model to both spaceborne and airborne datasets alludes to the possibility of sampling SSGM via Unmanned Aerial Vehicles (UAVs) (Pina and Vieira, 2022). Imaging instruments can be coupled with UAVs to achieve greater sampling frequency at a finer spatial resolution, while also being less susceptible to limitations of cloud cover, thus greatly complementing existing remote sensing platforms (Li et al., 2023). The potential of leveraging commercial UAV products to conduct sea surface imaging has been demonstrated by Wójcik-Długoborska et al. (2022) to study glacial suspended sediment plumes in the WAP.

Although this study demonstrates a new method for quantifying SSGM and its applicability across multiple remote sensing platforms, we stress that this is the first data product for quantifying SSGM remotely. While the SSGM predictions are validated with one of the most comprehensive in-situ δ¹⁸O datasets compiled for the WAP, there were no concurrent in-situ optical measurements; accordingly, our method takes an applied approach to correlate ocean color signals with field SSGM measurements. For the development of a first-generation model, we chose this applied approach for rapid implementation, so the results can be utilized by a wide range of end users for future projects. Additional laboratory studies are needed to understand the inherent and apparent optical properties of pure and diluted glacial meltwater.

The ability to monitor SSGM is directly relevant to understanding ecological dynamics that are being impacted by accelerating change along the WAP. The presence of SSGM is also found to coincide in regions of high chl-a concentrations. For example, phytoplankton abundance is found to be significantly correlated with glacial meltwater over the WAP shelf. This relationship extends >100km offshore and is persistent across years (Dierssen et al., 2002). The same dynamics between glacial meltwater and phytoplankton were also observed in a WAP fjord (Pan et al., 2020). In this study, we found similar spatial distributions between chl-a concentrations and SSGM but with clear distinctions between the two ( Figure 7 ), consistent with prior studies. These results have important implications for polar ecosystem research. Remote sensing studies of the WAP ecosystem have primarily relied on measurements of chl-a, SST, sea ice properties, and salinity (IOCCG, 2015). However, in recent years, glacial meltwater has been identified as an additional environmental variable that is important to the ecosystem in the WAP and broader polar regions (Dierssen et al., 2002; Arrigo et al., 2017; Meire et al., 2017; Pan et al., 2020). Freshwater addition experiments in Potter Cove found that a gradual decrease in salinity can shift phytoplankton communities dominated by large centric diatoms to small pennate diatoms, suggesting a phytoplankton response to low salinity is species-specific (Hernando et al., 2015). Similar results have also been found in the field over the WAP shelf where cryptophyte abundance was observed to coincide with relatively low salinity waters (S ≤ ~33.6 PSU) (Moline et al., 2004). In addition, other studies found that glacial meltwater’s impact on phytoplankton communities is likely beyond changes in salinity. In WAP fjords, glacial meltwater is also a pathway for delivering macro- and micro-nutrients to the surface. Sporadic but prolonged katabatic wind events bring up deep nitrate-rich water which is propagated along the surface with SSGM from the ice-ocean interface (Ekern, 2017; Pan et al., 2020). Dissolved iron is supplied to the surface via a similar mechanism, but it can also directly enter the ambient water column via sub-glacial and sub-marine melting (Sherrell et al., 2018; Forsch et al., 2021). These results indicate that glacial meltwater plays an important role near the ice-ocean boundary, and SSGM export from fjords likely serves as an important nutrient source for the broader WAP ecosystem (Forsch et al., 2021). In short, detection of SSGM from remote sensing can provide a new proxy for monitoring the physical environment, ecosystem, and their variabilities in the WAP.