3.1 Sample collection and preparation

During the 2017–2018 austral summer, 220 surface soil samples (∼ top 5 cm) were collected from 11 distinct ice-free areas (Roberts Massif, Schroeder Hill, Mt. Augustana, Bennett Platform, Mt. Heekin, Thanksgiving Valley, Taylor Nunatak, Mt. Franke, Mt. Wasko, Nilsen Peak, and Mt. Speed) along the Shackleton Glacier, including a subset of 27 samples previously analyzed for S, N, and O isotopes in nitrate and sulfate (Diaz et al., 2020b). At each area, we collected samples in transects (ranging from ∼ 200 to ∼ 2000 m in length) to maximize the geochemical variability. Our transects were also designed to capture the LGM transition, with some soils exposed throughout the LGM and others exposed following glacier retreat. GPS coordinates and elevations were recorded with each sample and later used to estimate the distance from the coast and distance from the glacier (defined as linear distance from the nearest glacier – Shackleton, tributary, or alpine). Once collected, the samples were stored and shipped frozen (−20^∘) to the Ohio State University.

Prior to geochemical analysis, the samples were dried at 50^∘C for at least 72 h, with the loss in mass attributed to soil moisture content. The dried soils were leached at a 1:5 soil to deionized water (DI) ratio, and the leachate was filtered through 0.4 µm Nucleopore membrane filters (Diaz et al., 2018, 2020b; Nkem et al., 2006). Due to the low sediment to water ratio, this leaching technique only dissolves the more water-soluble salts (Toner et al., 2013). These include salts with ClO${}_{\mathrm{4}}^{-}$, NO${}_{\mathrm{3}}^{-}$, Cl⁻, SO${}_{\mathrm{4}}^{\mathrm{2}-}$, ClO${}_{\mathrm{3}}^{-}$, and CO${}_{\mathrm{3}}^{\mathrm{2}-}$ + HCO${}_{\mathrm{3}}^{-}$. Process blanks were generated and analyzed to account for any contamination from the leaching process.

3.2 Analytical analysis of water-soluble anions, cations, and nutrients

The analytical techniques used here are similar to those reported by Diaz et al. (2020b). In brief, the analytes included anions (F⁻, Cl⁻, Br⁻, and SO${}_{\mathrm{4}}^{\mathrm{2}-}$) which were measured on a Dionex ICS-2100 ion chromatograph, cations (K⁺, Na⁺, Ca²⁺, Mg²⁺, and Sr²⁺) which were measured on a PerkinElmer Optima 8300 inductively coupled plasma optical emission spectrometer (ICP-OES), and nutrients (NO${}_{\mathrm{3}}^{-}+$ NO${}_{\mathrm{2}}^{-}$, PO${}_{\mathrm{4}}^{\mathrm{3}-}$, H₄SiO₄, and NH₃) which were measured on a Skalar San$++$ automated wet chemistry analyzer at the Ohio State University. Perchlorate (ClO${}_{\mathrm{4}}^{-}$) and chlorate (ClO${}_{\mathrm{3}}^{-}$) were measured using an ion chromatograph–tandem mass spectrometry technique (IC-MS/MS) at Texas Tech University (Jackson et al., 2012, 2015). All analytes are reported as listed. Total water-soluble salt concentration was calculated as the sum of all measured cations and anions. The precision of replicated check standards and samples was typically better than 10 % for all major anions, cations, and nutrients and better than 20 % for perchlorate and chlorate. Accuracy was typically better than 5 % for all major anions, cations, and nutrients, as determined by the NIST 1643e external reference standard and the 2015 United States Geological Survey (USGS) interlaboratory calibration standard (M-216), and better than 10 % for perchlorate and chlorate, as determined by spike recoveries. Precision and accuracy for individual analytes are located in Table S1. Detection limits for the analytes have been previous reported (Diaz et al., 2018; Jackson et al., 2012).

3.3 Data interpolation and machine learning

Inverse distance weighted (IDW) interpolations were performed for Bennett Platform, Thanksgiving Valley, and Roberts Massif using the Geostatistical Analyst tool in ArcMap 10.3. Since IDW is a deterministic method where unknown values are predicted based on proximity to known values, we chose those three sites as they had the most defined transects and relatively higher sample density. The interpolation parameters were constant with a power of 4, maximum neighbors of 15, minimum neighbors of 5, 4 sectors, and a variable search radius. These parameters were chosen such that they optimize for the lowest mean absolute error (MAE).

Multiple linear regressions were generated for all geochemical analytes except H₄SiO₄ (total of 15 dependent variables) with latitude, longitude, elevation, distance from the coast, distance from the glacier, and soil moisture as independent variables using built-in functions in R 3.6.3 (R Core Team, 2020). Random forest regression models were similarly generated using the RandomForest library. The random forest model is a machine learning algorithm that utilizes supervised learning algorithms to predict values given input predictor variables (Breiman, 2001). Multiple decision trees are run in parallel with a randomized subset of predictor variables, and the aggregate result of each tree is used to generate a predicted outcome. Since each tree is generated using a random sample and random predictor variables, the random forest model is effective at minimizing overfitting and handling outliers (Breiman, 2001). For both models, all geochemical data were log-transformed to ensure the data were normally distributed (verified using a Jarque–Bera normality test). Missing values were input as NA (not available).

Machine learning algorithms are widely used in a variety of disciplines from finance (Patel et al., 2015) to ecology (Davidson et al., 2009; Peters et al., 2007; Prasad et al., 2006) for both data prediction (regression) and classification. Recently, these techniques have been used for earth science applications, including geologic mapping (Heung et al., 2014; Kirkwood et al., 2016), air quality monitoring (Stafoggia et al., 2019), and water contaminant tracing (Tesoriero et al., 2017). We developed a novel application of machine learning to predict concentrations and gradients of water-soluble salts in Antarctic soils given set geographic parameters, similar to the approaches developed for stock market and real estate predictions (Antipov and Pokryshevskaya, 2012; Patel et al., 2015).

For our random forest models, any sparse missing values in Table S2 were estimated by averaging the geochemistry of the samples collected immediately before and after in the same transect. The new imputed dataset was split into a training set representing 86 % of the data (n = 189; Table S3) and a testing set representing the remaining 14 % (n = 31; Table S4) based on ideal model parameters described by Breiman (2001). The training dataset was used to generate the random forest models for each analyte. Each of the models was run with 2000 decision trees (ntree = 2000) to minimize the mean squared error. The number of random variables used for each node split in the decision trees was set to the recommended regression default of one third of the total number of variables to optimize the model randomness, which in our case was 2 (mtry = 2), following parameters described previously (Breiman, 2001). The scripts developed for both the multiple linear regression and random forest models are included in the Supplement.

4.1 Geochemistry of upper, middle, and lower zones

The maximum, minimum, mean, standard deviation, and coefficient of variation are reported in Table 1 for the measured geographic and geochemical data. Concentrations of water-soluble ions span up to 5 orders of magnitude and are variable across the region. Elevation, distance from the coast, distance from the glacier, and soil moisture are also variable and span up to 3 orders of magnitude. The highest elevation samples (>2000 m a.s.l.) were collected from Schroeder Hill, and the greatest soil moisture content is from Mt. Wasko at 12.3 % with a mean of 2.1 % for all samples.

Table 1Overview of geography, soil moisture, and water-soluble ions from the Shackleton Glacier region. The minimum values reported are those within the detection limits. Individual sample concentrations are detailed in Table S2.

SD: standard deviation, CV: coefficient of variation.

Download Print Version | Download XLSX

Shackleton Glacier region surface soils can be separated into three zones based on their water-soluble geochemistry: an upper zone near the Polar Plateau, a middle zone near the center of the glacier, and a lower zone where the glacier flows into the Ross Sea (Figs. 1, 2). The upper zone samples are characterized by the highest total water-soluble salt concentrations, with the highest values greater than 80 000 µg g⁻¹ at Schroeder Hill, while the lower zone samples have the lowest total salt concentrations, with the lowest values near 10 µg g⁻¹ at Mt. Wasko (Fig. 2a–c). The middle zone has intermediate values. Water-soluble N:P molar ratios generally follow a similar trend (Fig. 2d–f). The lowest N:P ratios are in the lower zone soils, while the middle and upper zones have more variable values. Concentrations of ClO${}_{\mathrm{4}}^{-}$ and ClO${}_{\mathrm{3}}^{-}$ follow similar trends as the total salts with less distinction between the middle and upper zones, though most concentrations in the lower zone are below the detection limit (Fig. 2g–l; Table S2).

Figure 2Total water-soluble salts, water-soluble N:P molar ratio, and ClO${}_{\mathrm{4}}^{-}$ and ClO${}_{\mathrm{3}}^{-}$ concentrations (log scale) were compared to elevation, distance from the coast, and distance from the glacier for samples from the three geographic zones (blue for upper, yellow for middle, gray for lower zones). Linear regression lines are plotted, for which dashed lines represent regressions where p>0.05 with a Bonferroni correction, and R² values are reported for each relationship. The horizontal orange lines represent nematode salt tolerance of ∼ 2600 (Nkem et al., 2006), and the green lines represent the Redfield ratio, N:P = 16, for phytoplankton in the ocean.

Download

Observed trends between the zones appear to be driven, at least partially, by geography. Regressions of total water-soluble salt concentration, water-soluble N:P ratio, and ClO${}_{\mathrm{3}}^{-}$ concentration with elevation, distance from the coast, and distance from the glacier are all positive (Fig. 2). The strongest relationships are between total salts and elevation, and N:P ratio and elevation with R² values of 0.59 and 0.52, respectively, and p values <0.001 with a Bonferroni correction which was applied to minimize the Familywise Type 1 error rate associated with multiple comparisons (Fig. 2a, d). The weakest relationships are between ClO${}_{\mathrm{4}}^{-}$ and distance from the coast and ClO${}_{\mathrm{3}}^{-}$ and distance from the glacier with R² values of 0.11 and 0.06, respectively (Fig. 2h, i). Distance from the glacier varies widely between individual zones with frequent overlaps, but there appears to be a moderate relationship with N:P ratio and total salts (Fig. 2c, f). Overall, total salt concentration has the strongest relationship with geography and ClO${}_{\mathrm{4}}^{-}$ has the weakest relationships.

Ternary diagrams highlight the specific geochemical gradients within and between the zones (Fig. 3). The anion ternary diagram only includes SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, and Cl⁻, which are the major water-soluble salts in the region (Claridge and Campbell, 1968; Diaz et al., 2020b). Though carbonate and bicarbonate salts have been identified in both lacustrine sediments and soils in Antarctica, previously measured concentrations in the Shackleton Glacier region were low, ranging from 0.07 % to 2.5 %, and bicarbonate salts were not identified in the highest elevation and furthest inland soils (Claridge and Campbell, 1968; Diaz et al., 2020b; Lyons et al., 2016). The most abundant anion for the upper zone is SO${}_{\mathrm{4}}^{\mathrm{2}-}$, which is greater than 99 % of the total anion budget in some Schroeder Hill and Roberts Massif samples, though other locations are dominated by NO${}_{\mathrm{3}}^{-}$ (Fig. 3). The anions are more evenly distributed in the middle zone, though the majority of samples are most abundant in NO${}_{\mathrm{3}}^{-}$ and Cl⁻. The lower zone has much lower SO${}_{\mathrm{4}}^{\mathrm{2}-}$ fractions than the upper zone, and the dominant anion is generally Cl⁻. The cation distribution is very similar for all three zones (Fig. 3). Na${}^{+}+$ K⁺ is the most abundant cation pair representing over 90 % of the total cations for many upper and middle zone samples, while Ca²⁺ is the second most abundant. In general, Mg²⁺ is the least abundant cation across all sampling locations.

Figure 3Anion and cation ternary diagrams for the three geographic zones.

Download

4.2 Statistical geochemical variability

A principal component analysis (PCA) using the correlation matrix (i.e., scale = TRUE) was performed in R (using factoextra; Kassambara and Mundt, 2017) and built in R software libraries to determine which geochemical variables most strongly differ across the samples. For the PCA, the first two principal components (PC1 and PC2) account for over 50 % of the total dataset variability at 44.2 % and 11.6 %, respectively. The different zones are correlated with different principal components (Fig. 4). The samples from the middle zone are positively correlated with PC1 and PC2. In the biplot, they plot in the upper-right quadrant with high concentrations of Cl⁻, NO${}_{\mathrm{3}}^{-}$, water-soluble N:P ratio, and Ca²⁺ and with a minor influence from soil moisture and H₄SiO₄. The upper zone samples generally plot along PC1 and are most associated with Sr²⁺, SO${}_{\mathrm{4}}^{\mathrm{2}-}$, Mg²⁺, Na⁺, K⁺, F⁻, ClO${}_{\mathrm{4}}^{-}$, and ClO${}_{\mathrm{3}}^{-}$. The samples from the lower, more coastal zone are negatively correlated with PC1 and are distinguished by their higher PO${}_{\mathrm{4}}^{\mathrm{3}-}$ concentrations. Most samples from all locations plot within the 95 % confidence interval ellipses. However, there are two strong outliers from Schroeder Hill and Mt. Heekin.

Figure 4Principal component analysis (PCA) biplot generated in R using factoextra and built in R software libraries with all anions, cations, nutrients, and soil moisture for the three geographic zones. The PCA is based on the correlation matrix (i.e., scale = TRUE). Principal component 1 and principal component 2 are plotted on the x and y axes, respectively. Shaded ellipses represent 95 % confidence intervals.

Download

Similar to the PCA, we performed a simple Spearman's rank correlation for the entire dataset to visualize the statistical dependence between all variables. Since a goal of this study is to relate water-soluble ion concentrations to geography, we focused on latitude, longitude, distance from the coast, distance from the glacier, and soil moisture. The strongest correlation coefficients are between Cl⁻ and latitude, elevation, and distance from the coast and Sr²⁺ and soil moisture (Fig. 5). Most other correlations are moderate to weak, though there appear to be notably stronger correlations between ClO${}_{\mathrm{3}}^{-}$ and latitude and distance from the coast, Ca²⁺ and longitude, elevation, and distance from the coast, NO${}_{\mathrm{3}}^{-}$ and latitude, and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and distance from the glacier. Longitude, elevation, and distance from the coast have the greatest number of strong and moderate correlations with the geochemistry data. Outside of the geographic parameters, Na⁺ is highly correlated with total water-soluble salts, likely representative of the high Na${}^{+}+$ K⁺ percentages (Fig. 3), and Sr²⁺ is highly correlated with K⁺, likely reflecting a common ion source.

4.3 Spatial interpolation and machine learning model performance

The total salt concentrations of individual samples at Bennett Platform produce the most defined interpolation gradient from the glacier front to further inland compared to Roberts Massif and Thanksgiving Valley (Fig. 6). Bennett Platform also has the smoothest salt concentration contours suggesting that the interpolation model is the strongest and most robust at this location. The second strongest interpolation is Thanksgiving Valley. Contrary to the measurements at Bennett Platform, Thanksgiving Valley has the highest salt concentrations in the center of the valley with lower concentrations to both the east and west. The lowest concentration contours are closest to the glacier for both Bennett Platform and Thanksgiving Valley, which is likely related to glacial history since the soils near the glacier are relatively younger than those further inland based on meteoric ¹⁰Be data (Diaz et al., 2020a). The interpolation from Roberts Massif does not have a distinguishable spatial trend.

The multiple linear regression and random forest models vary in their strength for the individual analytes. The highest R² value from the linear regression is 0.65 for Na⁺, while total water-soluble salts, water-soluble N:P ratio, ClO${}_{\mathrm{4}}^{-}$, and ClO${}_{\mathrm{3}}^{-}$ have values of 0.61, 0.60, 0.44, and 0.55, respectively (Table 2). The lowest R² value is for PO${}_{\mathrm{4}}^{\mathrm{3}-}$ at 0.17. The p values for all analytes are ≪0.001 even with a Bonferroni correction. The highest out-of-the-box explained variance values from the random forest models are for total salts and ClO${}_{\mathrm{3}}^{-}$ at 75 % and 63 %, respectively. The lowest explained variance is for Sr²⁺ at 37 %. The values of N:P ratio and ClO${}_{\mathrm{4}}^{-}$ are 52 % and 48 %, respectively. We also evaluated the most important and least important variables in the random forest models based on node purity. The most important variable for the majority of analytes is elevation, while distance from the glacier is most important for N:P ratio and latitude for ClO${}_{\mathrm{3}}^{-}$ (Table 2). The least important variables are distance from the coast and latitude for every analyte except ClO${}_{\mathrm{3}}^{-}$, for which distance from the glacier is least important.

Table 2Out-of-the-box multiple linear regression and random forest model statistics generated in R. All geochemical data were log-transformed.

Download Print Version | Download XLSX

5.1 Implications for ecological habitat suitability

By establishing geochemical zones for the Shackleton Glacier region, we can better understand the relationship between geochemistry and geography and ultimately biogeography. As stated in the introduction, we focused particularly on total water-soluble salt concentrations, water-soluble N:P ratios, and ClO${}_{\mathrm{4}}^{-}$ and ClO${}_{\mathrm{3}}^{-}$ concentrations.

5.1.1 Elevation and moisture controls on total water-soluble salt gradients

The elevational trends of total salt concentrations at the Shackleton Glacier are similar to those previously described in the TAM, where soils from higher elevation sites typically have higher salt concentrations (Bottos et al., 2020; Lyons et al., 2016; Magalhães et al., 2012). Our results are also consistent with those from Scarrow et al. (2014) who found that salt concentrations typically decreased with distance from the glacier in the Beardmore and Lennox-King glacial regions. Our total water-soluble salt interpolation maps highlight the spatial variability in Shackleton Glacier region soils (Fig. 6). The most spatially variable location is Robert Massif, which does not appear to follow local elevational, latitudinal, and/or distance inland gradients. This heterogeneity is not necessarily due to currently active soil leaching as the soil moisture values are not drastically different between the samples (Table S2). Though the variability in cation concentrations is likely due to the weathering of tills, scree, and bedrock (Claridge and Campbell, 1968), recent work on the isotopic composition of water-soluble nitrate and sulfate, the major anions in the upper zone, suggests a common, atmospheric source (Diaz et al., 2020b).

Figure 5Spearman's rank correlation matrix generated in r using the corrplot library. The colors represent correlation coefficients, indicating the strength and magnitude of the correlation. The blue box indicates the geographic variables and soil moisture, which were variables used in the multiple linear regression and random forest models. Familywise Type 1 error corrections were not applied in this analysis.

Download

Figure 6Inverse distance weighted (IDW) interpolations of total salt concentration for Roberts Massif (a), Bennett Platform (b), and Thanksgiving Valley (c). The color scale represents the 10 natural breaks in the data. Interpolations were created and mapped using the Geostatistical Analyst tool in ArcMap 10.3.

We argue that the heterogeneity in the total salt concentrations at Roberts Massif (Figs. 2, 6) is probably related to different and complex wetting history, in which seasonal snow patch melt may pool in local depressions, transporting water-soluble salts from slightly higher elevations and/or from saline wet patches (Levy et al., 2012). This is demonstrated on a larger scale at Thanksgiving Valley, a glacially carved valley where the higher concentrations of salts in the center of the valley are likely due to the transport of salts from nearby higher elevation slopes during melting events. This is further evidenced by the presence of two small, closed-basin ponds in the center of the valley which likely formed from glacial melt and may have been larger in size in the recent past (Diaz et al., 2020). Similarly, streams and meltwater tracks in the MDV leach soils and carry salts into closed-basin, brackish to hyper-saline lakes where salts are cryoconcentrated over time (Lyons et al., 1998). Our results suggest that elevation and wetting history are important contributors to total salt gradients in the Shackleton Glacier region as they influence the accumulation of salts and subsequent leaching from soils.

5.1.3 High and variable ClO${}_{\mathrm{4}}^{-}$ and ClO${}_{\mathrm{3}}^{-}$ concentrations

Our ClO${}_{\mathrm{4}}^{-}$ and ClO${}_{\mathrm{3}}^{-}$ concentrations include the highest measured in Antarctica to date and are comparable to concentrations from the Atacama and Mojave deserts (Jackson et al., 2015). Though not a strong correlation, the highest elevation samples (upper zone) have the highest ClO${}_{\mathrm{4}}^{-}$ and ClO${}_{\mathrm{3}}^{-}$ concentrations (Fig. 2g, j). Similar to NO${}_{\mathrm{3}}^{-}$, ClO${}_{\mathrm{4}}^{-}$ and ClO${}_{\mathrm{3}}^{-}$ are derived from atmospheric deposition, and, because of their high solubilities, their accumulations are related to wetting and glacial histories (Jackson et al., 2016, 2015). Therefore, soils which have been exposed for long periods of time and have not experienced snow or ice melt, such as those from Schroeder Hill and Roberts Massif, are able to accumulate high concentrations of ClO${}_{\mathrm{4}}^{-}$ and ClO${}_{\mathrm{3}}^{-}$. Interestingly, our ClO${}_{\mathrm{4}}^{-}$ concentrations are lower – maximum of ∼ 1.9 g L⁻¹ – than the highest recorded tolerance – 1.1 M (∼ 130 g L⁻¹) NaClO₄ – for the extremotolerant bacteria Planococcus halocryophilus, yet a recent study shows no detectable biomass for Schroeder Hill samples (Dragone et al., 2020). (Per)chlorates are strong oxidizers and are well established as toxic; thus the concentrations of ClO${}_{\mathrm{4}}^{-}$ and ClO${}_{\mathrm{3}}^{-}$ might be additional, crucial indicators of habitat suitability. However, the concentrations are highly heterogenous across our sampled locations (Fig. 2k–l), and unlike ClO${}_{\mathrm{3}}^{-}$, neither the multiple linear regression nor random forest models were able to adequately capture the variability in ClO${}_{\mathrm{4}}^{-}$ concentrations (Table 2).

Table 3Multiple linear regression and random forest statistics between predicted and measured concentrations from the Shackleton and Darwin Glacier regions. R² and p values are reported for the correlations between measured and predicted concentrations. Regression slopes and test statistic values (t) were calculated using the smatr library (Warton et al., 2012) in R to evaluate the null hypothesis (H₀) of slope = 1. Higher test statistic values (closer to 1) indicate that we reject the null hypothesis. All geochemical data were log-transformed.

^* t<0.5; ${}^{**}$ t<0.20.

Download Print Version | Download XLSX

Table 4Multiple linear regression and random forest mean absolute error (MAE) and root mean squared error (RMSE). All geochemical data were log-transformed for the analysis.

Download Print Version | Download XLSX

5.2 Machine learning as a tool to predict soil geochemical trends

We sought to evaluate our multiple linear regression and random forest models using a testing dataset from the Shackleton Glacier region (n = 31) and a second dataset from the Darwin Mountains (∼ 80^∘ S) (n = 10) (Magalhães et al., 2012). Few published/available TAM datasets include sample GPS coordinates, soil moisture, and water-soluble ion geochemistry. As stated in Sect. 3.3, the Shackleton Glacier region test data were not included in the random forest model generation so we could evaluate our models with an independent dataset. For the Darwin dataset, distance from the glacier, distance from the coast, and elevation were determined using the Reference Elevation Model of Antarctica (REMA), while location, soil moisture, and geochemistry were retrieved from the literature (Howat et al., 2019; Magalhães et al., 2012). We evaluated all 15 analytes from the original models with the Shackleton dataset and, due to a lack of data, only evaluated seven analytes from the Darwin soils (Fig. 7).

Figure 7R² values for the multiple linear regression and random forest predicted and measured values for the different analytes (Table 3). Test datasets include the Shackleton Glacier region (n=31) and the Darwin Glacier region (n = 10) (Magalhães et al., 2012). Analytes with slopes near 1, indicating good agreement between measured and predicted values, are indicated (^* t<0.5; ${}^{**}$ t<0.20).

Download

Both the multiple linear regression and random forest model outputs are moderately well-correlated for the Shackleton dataset, as determined by Pearson correlations between the measured and predicted values (Fig. 7a; Table 3). The random forest models outperform the linear regression models for nearly every analyte with the exception of Sr²⁺, NH₃, and PO${}_{\mathrm{4}}^{\mathrm{3}-}$, and nearly all p values are <0.001. For Cl⁻, in particular, the random forest model significantly outperforms the multiple linear regression model with R² values of 0.67 and 0.16, respectively. N:P molar ratio is the most accurately predicted analyte with R² values of 0.88 and 0.59 for the random forest and linear regression models, respectively. However, the highest R² value for the multiple linear regression model is for Na⁺ at 0.64 (Table 3). In terms of our analytes of interest regarding habitat suitability, total salts have the second strongest correlation (following N:P ratio) with the random forest model (R²=0.81), followed by ClO${}_{\mathrm{4}}^{-}$ (R²=0.78) and ClO${}_{\mathrm{3}}^{-}$ (R²=0.74). MAE and root mean squared error (RMSE) values indicate that the random forest models also have a smaller error compared to the multiple linear regression models (Table 4). MAE values are lower than RMSE values for both models, indicating the strong influence of outliers in the testing dataset. This is unsurprising as the standard deviation and coefficient of variation values for the entire dataset are relatively large for all analytes. Additionally, the outliers are likely one reason why the random forest models are stronger than the multiple linear regression models.

Similar to the model performance in the Shackleton Glacier region, the water-soluble ion predictions for the Darwin Glacier region are more strongly correlated with measured values in the random forest models compared to the multiple linear regressions (Fig. 7b). In fact, the linear regression models fail for nearly all the Darwin samples, and most concentration outputs are negative, which is likely due to overfitting during model generation. Here, Ca²⁺ and K⁺ are exceptions, and the multiple linear regression models outperform the random forest models in both cases. MAE and RMSE values for both models are higher than those for the Shackleton dataset (Table 4). On the other hand, the random forest models perform particularly well for some analytes. Though a small sample size, the R² values for N:P molar ratio and Ca²⁺ are 0.68 and 0.66, respectively, with p values ≪0.001. Total salts is moderately correlated (R²=0.47). It is unclear why some analytes, such as N:P molar ratio, are the most accurately predicted, though we suspect that this is due to (1) weathering trends of local lithology across the TAM since chemical weathering is probably the major source of these ions and (2) deposition and accumulation of atmospherically derived ions at higher elevations (Diaz et al., 2020b).

It should be noted that the R² values simply measure the strength of the correlations between the measured and predicted values. We performed slope tests by fitting bivariate lines using the standardized major axis (SMA) to further understand the relationship between the two values using the smatr library in R (Warton et al., 2012). For this test, we specifically evaluated the null hypothesis (H₀) in which slope = 1, which would indicate whether an ideal, direct 1:1 relationship exists between the measured and predicted values. Test statistic values (t) were used to measure the sample correlation between the residuals and fitted values (Warton et al., 2012). Test statistic values near 1 indicate that we reject the null hypothesis. In other words, higher absolute test statistic values indicate a slope other than 1. Of the 15 analytes in the Shackleton dataset, 5 analytes have slopes near 1 for the multiple linear regression models and 11 for the random forest, as indicated by test statistic values less than 0.5. For the Darwin, no analytes have test statistic values less than 0.5 (Fig. 7; Table 3).

These data indicate that while some analytes have high correlations between measured and predicted values, the models perform best with the Shackleton Glacier region soils. Additionally, though the relationship may not be 1:1, the random forest models are effective at predicting the measured geochemical gradients. For example, similar to our data, the Darwin Glacier samples generally have greater water-soluble N:P ratios and total water-soluble salt concentrations further from the glacier and at higher elevations (Magalhães et al., 2012), a trend that is reflected by our model results despite offset values. Additionally, corrections for the offset of the model from a slope = 1 (i.e., multiplying the model output value by the regression slope) can be made to better estimate specific concentrations, though the difference between modeled and measured values can still be up to 2 times greater. Our sample size for building the multiple linear regression and random forest models is small. We anticipate that as more data are collected throughout the CTAM, these data can be added to the model training dataset, expanding our prediction capabilities and increasing model reliability.