Spatial modelling of available aluminum in glacial till: informing a boreal forest soil carbon predictor

Boreal forest soils store 30 % of global forest soil carbon, making them a crucial component of the carbon cycle. However, climate change at high latitudes is increasing temperatures and producing unpredictable precipitation patterns. Soil organic carbon (SOC) formation and stabilization is tied to precipitation patterns, thus climate change will inevitably influence the stability and longevity of boreal forest SOC. The current size and distribution of this reservoir is poorly understood, creating uncertainty under current and future climate scenarios. Previous research demonstrates mineral soil properties may be used to model boreal forest SOC accurately. The surface slope, depth of carbon enriched horizon, and climate characteristics are important parameters for modelling SOC and can generally be obtained or estimated via remote sensing. However, information about the availability of reactive aluminum – non-crystalline forms capable of interacting with organic matter to form stable carbon rich organometal complexes in mineral soils – is not widely available but is controlled by the weatherable aluminum in soil parent material. To bridge this gap, the Newfoundland and Labrador till geochemistry dataset was used here to map and model weatherable aluminum availability in glacial till across the island of Newfoundland as a function of geology and climate. The Random Forest Algorithm was employed to develop two models: one relying solely on the strength of geological and climatic variables, and the other drawing additionally on the spatial context of sample points. The first model performed well (R²=0.60), however, adding a spatial component increased the performance of the second model (R²=0.71). Bedrock type and proximity of the samples to certain units were indicated to be the strongest controls on aluminum availability, while environmental factors were less influential. Additionally, model uncertainty was calibrated empirically and mapped spatially, providing reliable and actionable information about the confidence of predictions over the study area. This project demonstrates the value of predictive geospatial modelling for till geochemistry mapping and delivers key aluminum availability predictions useful for future efforts in deriving SOC reservoir estimates across Newfoundland.

1 Introduction

Soils are a crucial component of the global carbon cycle, with soil carbon (C) stocks on average three times larger than the atmospheric C reservoir (Scharlemann et al., 2014). Of this, boreal forest soils store 30 % of global forest soil C alone (Scharlemann et al., 2014). This is primarily stored in podzols – an acidic mineral soil characterised by a nutrient- and organic-depleted eluvial (Ae) horizon overlying the enriched illuvial (B) horizon (Hagemann et al., 2010). Soil organic carbon (SOC) is sequestered when dissolved organic matter leached from the overlaying organic layer infiltrates into the mineral soil and interacts with reactive metals, precipitating as organometal complexes (OMCs) in an organic matter enriched soil horizon (Baldock and Skjemstad, 2000). Transport of SOC through the soil profile to the persistent deep mineral soil C pool depends on water infiltration, with SOC stabilisation influenced by climate regulated microbial composition and activity (Jackson et al., 2017; Wang et al., 2022). Soil C storage potential in boreal forests is linked to a variety of soil parameters such as texture, clay content, acidity, cation availability, and environmental factors such as climate, topography, and parent material (Jenny et al., 1949; Oades, 1988; Slessarev et al., 2022; Patrick, 2023). However, in the young podzolic soils that dominate most boreal forests, reactive metal availability, dominated by reactive Al availability (Patrick et al., 2022), are one of the strongest predictors of SOC (Rasmussen et al., 2018).

Amplified climate change in high latitudes raises concerns about vulnerable boreal SOC stocks as the current and future size and distribution of this reservoir under various climate scenarios are poorly constrained, introducing uncertainty in models relying upon such estimates (Crowther et al., 2016; Friedlingstein et al., 2014; Wang et al., 2022). Consequently, unravelling climate-SOC feedbacks is crucial to predict the SOC stock evolution under future climate scenarios (Friedlingstein et al., 2014; Jackson et al., 2017). Additionally, evaluating this reservoir at global, regional, and local scales is key to mitigating climate change. Land use and resource extraction management can mitigate terrestrial CO₂ loss into atmospheric or aquatic reservoirs (Wang, 2024). Such practices include modifying agricultural techniques to include low tillage practices and fertilizer amendments, or reforestation and harvesting plans that encourage SOC retention (Hillier et al., 2012; Ameray et al., 2021).

Friedlingstein et al. (2014) analysed global SOC estimates from 11 participants in the 5th Climate Model Intercomparison Project and found estimates ranged from 510 to greater than 3000 Pg, with wide error margins disproportionately affecting high latitude regions. This is because SOC is modelled globally using earth system models (ESMs), which synthesize results of smaller models to provide insight into environmental context, climatic conditions, and soil properties (Pierson et al., 2022). Podzol-specific SOC storage mechanisms are often not incorporated into such models, and boreal regions are thus poorly modelled – likely a consequence of low sample density in high latitude areas (Paré et al., 2024). For example, productivity attributes often dictate SOC predictions (Zhang et al., 2020), however these are unreliable at high latitudes where organic input is scant and undergoes protracted decomposition due to low temperatures and elevated acidity (Paré et al., 2024; Zhou et al., 2024). Additionally, in some cases soil texture properties have been improperly positively related to SOC storage; the CENTURY climate model (Parton et al., 1993) associates higher clay fractions with increased SOC stocks, obscuring the importance of other mineralogical controls (Rasmussen et al., 2018; Zhang et al., 2020; Patrick, 2023; Paré et al., 2024). Similarly, the scale at which SOC is estimated is consequential. Research suggests that a fine scale approach delivers higher accuracy and more realistic estimates over regional extents (Deluca and Boisvenue, 2012; Jackson et al., 2017; McNicol et al., 2019; Paré et al., 2024; Pierson et al., 2022). For example, constraining models to the biome scale may normalize the effects of climate and illuminate SOC mechanisms that may otherwise be concealed (Patrick, 2023).

A growing body of knowledge suggests that metal oxyhydroxides serve as a valuable SOC predictor in boreal and other humid forest soils. For example, Spinola et al. (2026) demonstrate the association between podzolic SOC distribution and Fe and Al oxyhydroxides In soil, these forms of Fe and Al found in reactive poorly crystalline minerals (PCMs) stabilize carbon through organometal interactions including co-precipitation to form organometal complexes (OMCs) and adsorption to OMCs.The prevalence of these reactive metal forms in soil is a function of their content in weatherable mineral forms in the soil parent material (Slessarev et al., 2022), typically, the underlying glacial till in boreal regions. For example, higher Fe and Al in mineral soils were observed in boreal forests growing on highly weatherable ultramafic-rich till with high Al-rich feldspar content in contrast to similar forests overlaying more highly weathered, Al poor, sedimentary-based till (Patrick et al., 2022). Locally, the C saturation of organometal complexes (OMCs), the proportion of SOC in OMCs, and soil reactive Al content, (forms of Al in the soil that will directly participate in OMC formation) are the primary predictors for soil C in inaccessible, poorly sampled upland boreal forest soils on the Island of Newfoundland, further highlighting Al availability as a key predictor of SOC (Patrick, 2023). Therefore, Al availability – Al with the potential to contribute to reactive Al in soil – is a function of both Al content and mineral weatherability in soil parent material. The chemical index of alteration (CIA) which compares the content of highly stationary Al oxides to more mobile Ca, Na, and K oxides (Nesbitt and Young, 1984) is a commonly accessible approach to evaluating the degree of weathering a mineral has undergone. Combining soil parent material Al content with CIA provides an estimate of Al availability that can serve as a proxy for a key geologic control on SOC storage (organometal interactions), as reactive metals prove to be a more reliable indicator of SOC storage capacity in podzols than clay and silt (Rasmussen et al., 2018; Patrick et al., 2022; Spinola et al., 2026).

Previous work illustrates that slope, depth of enriched horizon, and climate can be modelled to extrapolate the saturation of OMCs and proportion of SOC within humid boreal forest soil OMCs (Patrick, 2023). However, soil Al % is not so easily estimated and is currently obtained through intensive field work and laboratory analyses. Assessing the controls influencing Al OMC precipitation is key to predicting SOC in these soils (Patrick et al., 2022; Patrick, 2023). Bedrock geological maps are often used as a proxy for parent material. However, this ignores the capacity of different lithologies to supply reactive mineralogical material to support soil forming processes through differences in weathering rates. Geological maps also explain little about the enduring effects of parent material on soil composition (Lawley et al., 2008). Unsurprisingly, surficial geology has been found to be a consistently significant control of SOC at local to regional scales (Adhikari et al., 2020). On the Island of Newfoundland, and in many boreal regions, parent material is dominantly glacial till, not bedrock, rendering it especially relevant to incorporate this dimension into SOC modelling. First mapping the till composition and using this to support a future sampling design to capture a range of Al % in soils across a variety of slopes and climatic conditions would provide a means to explore this parent material–Al OMC connection (Patrick, 2023).

The aim of this paper is to create a high-resolution geospatial model predicting the distribution of Al-rich weatherable minerals in glacial till in Newfoundland. This is a critical step toward attaining accurate regional estimates of boreal forest soil SOC stocks across the Island of Newfoundland. We model Al % spatially as a function of covariates believed to influence Al geochemistry. Because till geochemistry is controlled by geological parent material, and because Al is a relatively immobile element, we hypothesise that bedrock geology will be the strongest predictor for Al %. Two models are trialled in this study: a highly specific spatial model using both environmental covariates and geographic location information, and a more general aspatial model relying solely on environmental covariates. We hypothesise that the spatial model may be useful in areas that have been extensively sampled and the aspatial version may have a higher capacity to extrapolate predictions to areas with fewer samples. Predictive maps of Al and calibrated spatial estimates of their uncertainty are produced as objectives of this project. Calibration of predictive uncertainty is crucial for delivering accurate estimates of model confidence that can be used directly in future SOC modelling. These objectives will contribute to the establishment of a continuous spatial SOC stock estimate for the Island of Newfoundland.

2 Methods

2.1 Till dataset

Publicly available till geochemistry data was sourced from the Geological Survey of Newfoundland and Labrador (GSNL; Newfoundland and Labrador Geological Survey, 2024). This dataset was developed from till surveys conducted by the GSNL between 1985 and 2023 and spans much of the island (see Fig. 1). The inconsistent sample distribution can be attributed to sampling primarily being conducted according to National Topographic System map sheets. Samples were collected 40–60 cm deep on average and were large enough to obtain ∼0.50 kg of till matrix, with larger clasts discarded (Newfoundland and Labrador Geological Survey, 2024). Samples were dried and sieved using either a <180 or <63 µm sieve. Samples underwent gravimetric analysis and the remaining solid was digested in 15 mL of concentrated HF, 5 mL of concentrated HCl, and 5 mL of 50 volume percent HClO₄. Further addition of concentrated 10 % HCl was added if dissolution was not achieved. Geochemical composition of the solution was determined by the GSNL using inductively coupled plasma emission spectroscopy (ICP-OES). Results are given in wt % for elements of interest which were Na, Ca, K, and Al.

Figure 1 Distribution of till geochemistry samples from the Geological Survey of Newfoundland and Labrador across Newfoundland.

2.2 Covariates

Covariates were selected to represent geological setting (surficial geology, regional geology), landscape characteristics (slope, ruggedness, elevation), and environmental context (climate region, distance from ocean). Additionally, a second model was developed that included layers representing the distance to each bedrock unit to provide spatial context. This ensures the model will consider the relationship of till geochemistry to the immediate bedrock unit, and within the greater geological setting (Fig. 2).

Figure 2 Project workflow for aluminum availability (Al %) models 1 and 2 and their respective standard deviations (SD).

2.3 Preprocessing

Till samples with Al, Na, Ca, and K values below the detection limit were removed, as were those with extremely low arbitrary values, which were identified as artefacts of the ICP OES analysis. Elemental composition was reported in terms of wt % oxide. Soil organic carbon is stored in Al organometal complexes which are formed from fresh, primary weathering products rather than extensively weathered forms of Al in the underlying glacial till (Slessarev et al., 2022). To assess the degree of mineral weathering in each till sample, the chemical index of alteration (CIA; Nesbitt and Young, 1984) was calculated following Eq. (1):

$$\begin{array}{}\text{(1)}& \mathrm{CIA}={\displaystyle \frac{{\mathrm{Al}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{3}}}{{\mathrm{Al}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{3}}+\mathrm{CaO}+{\mathrm{K}}_{\mathrm{2}}\mathrm{O}+\mathrm{NaO}}}\times \mathrm{100}.\end{array}$$

While concentrations of stable Al may be higher in a more weathered sample due to loss of more mobile elements, it is less likely to contribute to the formation of soil Al oxyhydroxides, which are depleted with extensive weathering (Slessarev et al., 2022). Therefore, the available Al (Al %), or the fresh weatherable Al available to support Al oxyhydroxides in the overlying soil, was then estimated using both the Al₂O₃ content and the CIA through Eq. (2):

$$\begin{array}{}\text{(2)}& \mathrm{Al}\mathrm{\%}={\displaystyle \frac{\mathrm{100}-\mathrm{CIA}}{\mathrm{100}}}\times {\mathrm{Al}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{3}}\times \mathrm{100},\end{array}$$

where Al % is expressed as a weight percent of oxides.

Al % is a novel method for estimating the Al available for Al phases that can participate in OMC formation but has intrinsic limitations resulting from the use of the CIA and the nature of the till geochemistry dataset which reports values in terms of bulk elemental composition. Unlike selective dissolution extractions or other methods of directly measuring the operational Al in soil, Al % is more reflective of a broader Al pool available in the parent material. While encompassing the weatherable Al pool available in the till, and the source of Al supporting soil OMC formation, the Al % is ascertained using bulk elemental composition of till. This cannot parse existing reactive Al from Al still stored as a stable mineral phase yet to supply fresh material through future weathering. As such, the till Al % expressed here represents a proxy useful for informing the potential reactive Al content in the overlying soil, rather than a direct measure of the Al pool operational in SOC organometal interactions. Further limitations associated with use of the CIA are discussed in Sect. 4.3.

A set of environmental covariate rasters were developed to model Al % spatially (Table 1). Al % is intrinsically related to weathering, which is influenced strongly by water infiltration and climate. Climate covariates dictating water availability (precipitation, snowfall, temperature, and evapotranspiration) were selected to represent this information. Water infiltration is decisively affected by topographic factors resulting in the inclusion of slope, ruggedness, and elevation. Distance to the ocean was also included; following glacial retreat, marine inundation affected coastal areas resulting in reworking and weathering of surficial deposits near the coast. This form of weathering was thought to potentially impact Al % in these areas. Distance to bedrock layers were additionally developed using Euclidean distance from the nearest bedrock unit, which is potentially relevant since glacial till samples are not derived exclusively from the bedrock unit upon which they lie, but also from pulverised rock from nearby units, with distant units affecting the composition less than the nearby units.

GSNL data layers were all converted to rasters with a resolution of 20 m in ArcGIS Pro across the full extent of the Island of Newfoundland. Categorical variables were coded numerically for modelling (bedrock geology unit 1–46; surficial geology units 1–9; climate zones 1–7). The “Near” tool in ArcGIS Pro was used to obtain the distance from the coastline to the centre of each raster cell, represented as a point layer, which was rasterized to match the resolution and extent of the other layers. Mean monthly snowfall, evapotranspiration, and air temperature were acquired from the FLDAS Noah Land Surface Model at 0.1° resolution for the period 1 January 1982 to 30 April 2024 (McNally and NASA/GSFC/HSL, 2018). The mean annual precipitation layer was derived from the WorldClim layer obtained from the ArcGIS Living Atlas at 3 km resolution (Fick and Hijmans, 2017). The data were masked to the study area, converted to vector format, and joined to one single 20 m grid-like polygon layer with an extent matching the raster layers. Values were interpolated in areas missing data using nearest neighbours within 10 km. This was done in preparation for use of the “Multivariate Clustering” tool in ArcGIS Pro, which was then used to identify seven climate zones using the mean annual precipitation and mean monthly snowfall, evapotranspiration, and air temperature values from the grid-like layer. The map of the climate zones was then converted to a raster. The Newfoundland Digital Elevation Model (NLDEM) (Government of Newfoundland and Labrador Fisheries, Forestry, and Agriculture, 2023) was used for elevation values. Although the NLDEM has a higher resolution, it comprises rounded integer values making the calculation of terrain attributes from this data source unsuitable. Therefore, slope and roughness were calculated from the Canada high resolution digital elevation model (Government of Canada, 2024) in R using the MultiscaleDTM package (Ilich et al., 2023).

Table 1 Characteristics and sources of the geospatial data layers used for model development. n/a indicates that these are not applicable to a given layer.

^a resampling method: nearest neighbor, ^b resampling method: bilinear

Additional data layers were produced representing spatial and geological context of points. The Euclidean Distance tool in ArcGIS Pro was used to create layers indicating the nearest distance to each of the GSNL Regional Bedrock Geology units. This yielded an additional 37 raster layers that were used only in Model 2 (the spatial version).

2.4 Modelling

Till geochemistry data and all rasters were loaded into RStudio using the terra package (Hijmans, 2025). The rasters were stacked, and the values were extracted at all till sample locations and added to the till geochemistry data frame. Incomplete observations were removed, producing a dataset of 17 328 samples. The Regional Surficial Geology variable was reclassified to “drift poor” in instances of a till sample falling on bedrock.

A Random Forest model was constructed to predict the Al % observed in the till samples using the environmental geospatial covariates. The Random Forest algorithm is effective at limiting overfitting, integrating both categorical and numerical data simultaneously, and parsing complex relationships between covariates and response variables, which are desirable characteristics for modelling SOC at a regional scale (Breiman, 2001; Heung et al., 2014). Random Forest is a machine learning approach that uses an ensemble of decision trees to make predictions. Decision trees are grown from randomly selected “bootstrap” samples of the dataset, which each comprise approximately two thirds of the total sample size as a function of sampling with replacement. Decision trees split data observations (i.e., Al %) into increasingly homogenous groups – at each split assigning a variable (m) best explaining the division (Breiman, 2001). This m is selected from a randomly generated portion (mtry) of all total variables (M). Predictions are made based on an average of the ensemble, meaning the variables which best explain the Al content are used to make a prediction for each tree, and the average of these forms the final prediction. The samples that were not drawn for each decision tree (i.e., the “out-of-bag” samples) are used to estimate the model accuracy. This built-in validation avoids splitting the dataset into “train” and “validation” partitions during model training. Additionally, 10-fold cross validation was performed to test the model predictions. Finally, Random Forest contains functionality to evaluate relationships between predictor and response variables (i.e., between Al % and environmental variables).

Two models were developed to predict values of Al % across the island of Newfoundland. Model 1 (the aspatial model) included seven covariates; elevation, distance from the ocean, climate zone, bedrock unit, surficial geology unit, slope, and ruggedness. Model 2 (the spatial model) used the same covariates as Model 1 but did not include the distance to ocean layer and replaced the bedrock geology variable with the 37 individual distance to bedrock geology variables enabling the model to leverage spatial trends in the geochemistry data. This addresses the varying provenance of these materials, as glacial till composition is not only influenced by the underlying geological unit, but also from neighbouring lithologies, because of glacial transport. It should be noted that although both results were mapped spatially, the titles “spatial” and “aspatial” model do not refer to this. Instead, these titles reflect the inclusion of information about sample location in the second model, and the lack of this information in the first model. The Random Forest hyper-parameters for Model 1 were M=7 (total variables), mtry = 4 (number of variables to draw at each split), and ntree = 1000 (number of trees to grow). Hyper-parameters for Model 2 were M=44, mtry = 13, and ntree = 2000. The R² and RMSE were calculated to evaluate each model via out-of-bag cross-validation. The resulting predictions and the predicted standard deviations δ of each ensemble model were output as raster layers to visualize spatially the distribution of Al % and model uncertainty, respectively.

2.5 Spatial model uncertainty

Following initial modelling and predictions the calibrated uncertainty of the models was derived and mapped spatially to convey confidence in Al predictions. Calibration refers to using the out-of-bag model error to scale δ to more realistic values. This is a necessary step to ensure that predicted δ reflects known model error, increasing confidence in these predictions. To accomplish calibration of the model uncertainty, the data were split into a test and training set comprising of 10 % and 90 % of the data respectively, with the training data used to generate calibration values and test data withheld to validate the optimized calibration values. Uncalibrated estimates of uncertainty were extracted from the out-of-bag standard deviation of the predictions of all trees of the Random Forest ensemble. To ascertain whether these estimates of uncertainty required calibration, graphical visual analysis of the training data was first performed (Palmer et al., 2022). The standard deviation (δ) of the out-of-bag predictions and the RMSE of residuals were calculated for the training data (90 % of data) for both the spatial and aspatial models. The withheld test data (10 % of data) was used to evaluate the calibration in both cases. If the model prediction δ accurately represents the uncertainty of predictions, the ratio of the residuals to the δ (the r statistic) should be normally distributed with a standard deviation of 1 (Palmer et al., 2022). The r statistic was calculated and plotted as a histogram for both models. While this emphasizes the relationship between δ and RMSE, it does not conclude whether larger RMSE values correlate with larger δ values (Palmer et al., 2022). To evaluate the linear relationship between model δ and cross-validated error, the predictions were divided into 15 equal frequency bins by δ, and the RMSE of residuals from each bin is calculated. The RMSE and δ are plotted against an identity function and checked for linearity. Preliminary results of both graphical analyses of the training data indicate that while the model uncalibrated standard deviation predictions δ_UC estimates are fair, they may be improved with calibration.

The Random Forest model prediction uncertainty was thus calibrated according to Palmer et al. (2022), which assumes that the error of predictions – or, residuals R(x) from the out-of-bag data D_OOB – are linearly related to δ_UC for any given point, or ${\mathit{\delta }}_{\mathrm{Cal}}=a{\mathit{\delta }}_{\mathrm{UC}}\left(x\right)+b$. Calibration values were obtained by optimizing Eq. (3) using the Nelder-Mead algorithm in R:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}a,b& ={\mathrm{argmin}}_{a^{\prime },b^{\prime }}\sum _{x,y\in D_{\mathrm{OOB}}}\ln \left(\mathrm{2}\mathit{\pi }\right)\\ & +\ln \left({\left(a^{\prime }{\mathit{\delta }}_{\mathrm{UC}}\left(x\right)+b^{\prime }\right)}^{\mathrm{2}}\right)+{\displaystyle \frac{R(x{)}^{\mathrm{2}}}{{\left(a^{\prime }{\mathit{\delta }}_{\mathrm{UC}}\left(x\right)+b^{\prime }\right)}^{\mathrm{2}}}}.\end{array}\end{array}$$

This minimizes the negative log likelihood that residuals R(x) are normally distributed with mean=0, and $\text{SD}=a{\mathit{\delta }}_{\mathrm{UC}}+b$ to find optimized values of a and b (slope and intercept) that are used to calibrate δ_UC (uncalibrated standard deviation) via linear model. To validate these results the δ_C of the test data were then calculated by scaling δ_UC by a and b using Eq. (4):

$$\begin{array}{}\text{(4)}& {\mathit{\delta }}_{\mathrm{Cal}}=a{\mathit{\delta }}_{\mathrm{UC}}+b.\end{array}$$

These were graphically and statistically validated using the same methods detailed above. Following this, the calibration factors were applied to the entire predicted δ data and mapped spatially.

3 Results

3.1 Mapped predictions and uncertainty

Both models predicted high values of Al % in several places throughout the island, notably the Mount Peyton (central northeast), Avalon (most easterly), and Deer Lake (south central west) areas (Fig. 3a, b). Highs also extend into the southwestern portion of the island, near Port Aux Basques. Low values were predicted in south central Newfoundland, over known weathered sedimentary packages. These patterns were predicted consistently by both models and are not surprising given the geological context. While both models predicted a high Al % near Gros Morne (central western coast), the spatial Model 2 predicted high values extending into the Great Northern Peninsula (northwest) and generally predicted higher values than the aspatial Model 1 in the southwestern area of the island.

Figure 3 Aluminum availability (Al %) predictions and corresponding calibrated uncertainty predictions given as Al % (SD) across the island of Newfoundland with aspatial (a, c) and spatial (b, d) models.

The prediction uncertainty δ was highest closer to the ocean and with increasing distance from samples – especially on the Great Northern Peninsula (Fig. 3c, d) where samples are sparse or non-existent (Fig. 1). While this pattern was observed in both outputs, it was more prominent in the spatial model. The aspatial model showed δ spread relatively evenly throughout the island, with lower values in areas close to samples and areas near the coast or distant from samples showing higher values. Unsurprisingly, bedrock units with few samples also showed a high δ.

3.2 Uncertainty calibration

Calibrated uncertainty of both models conforms closely to the identify function and the standard normal distribution. Scatterplots or the r statistic clarified that for both models, uncalibrated points lie beneath the identity function, indicating that both models overestimated δ. They also indicated a positive slope between calibrated and uncalibrated values, demonstrating positive correlation between the size of residuals and the size of δ estimates. Deviations from this pattern exist – though primarily in poorly sampled bins (Fig. 4e and f). In both cases, calibration improved the fit of data to both the normal distribution and the identity function. After calibration, normality and accuracy of the r statistic was improved, with standard deviations calibrated from 0.46 and 0.41 to 0.97 and 0.94, and the RMSE of the δ compared to the out-of-bag errors decreased from 0.317 and 0.264 to 0.207 and 0.160 for the aspatial and spatial models respectively – a decrease of approximately 34 % and 47 %.

Figure 4 Distribution of the R statistic (a, b), distribution of sample δ (c, d), and residuals vs. δ plots of both raw data and binned data (e, f) for aspatial (a, c, e) and spatial models (b, d, f).

The calibration was most reliable in areas of moderate uncertainty. Figure 4e and f show the highest uncertainty bin deviating significantly from the identity function, indicating that samples with a calibrated δ of greater than approximately 0.22 and 0.25 for the aspatial and spatial models respectively are underestimated. These points represent approximately 7 % of samples from both models. Application of the calibration factors to the δ rasters translated to 19 % and 40 % of δ being underpredicted across Newfoundland aspatially and spatially, respectively (Fig. 5). These underpredicted areas coincide entirely with areas of low confidence due to a paucity of samples. This is to be expected; calibrating uncertainty addresses systemic error within model predictions, not error associated with issues in sample distribution. Figure 4a, b, e, and f clearly demonstrate that overall, the uncertainty calibration corrects the consistent underprediction of δ by both models.

Figure 5 Map of δ calibration accuracy across Newfoundland for aspatial (a) and spatial (b) models.

3.3 Variable importance

Cross-validation demonstrated that Al % in glacial till could be modelled accurately using environmental variables, and with lower uncertainty with the addition of spatial context. The aspatial and spatial models attained R² values of 0.59 and 0.71, and RMSE values of 0.34 Al wt % and 0.29 Al wt % respectively. These results suggest that the spatial model explained >10 % more of the variance within sample data than the aspatial model, corresponding to a prediction uncertainty which was lower, on average, by ∼0.12 Al wt %.

The most influential variables for the aspatial model were bedrock units, distance from the ocean, elevation, and unit size (Table 2). Climate zone, slope, ruggedness, and surficial deposit type had the least control on Al %. Distance to the ocean was not expected to exert this much influence, and could represent spatial autocorrelation between samples, which motivated the second spatial model (see Methods section). The importance % score of the bedrock variables suggests that while environmental variables can be used to model Al %, the till samples relationship to bedrock and to one another are the most important prediction factors. Partial dependence plots for the spatial model were reviewed during the modelling process but due to the large number of bedrock predictors the interpretation of the plots was convoluted. This is largely because individual bedrock units are of great importance to some points (locally), but not to others (those at a distance from the unit). The value of these plots was diminished due to this, and as a result they were not included in the study.

Table 2 Aspatial and spatial model variables of importance given as a % of total variance explained.

4 Discussion

4.1 Utility of the approach

Results illustrate the value of geospatial environmental variables and spatial context for mapping and modelling available Al. The resulting map addresses an underlying need for incorporation of geological controls into future SOC modelling in boreal forests and will directly inform regional scale modelling efforts as recommended for Newfoundland (Patrick et al., 2022; Patrick, 2023). This approach has utility in areas where existing soil samples are sparse, but geological information is abundant. Till sampling is widespread due to its applications in mineral exploration, and Al, Ca, Na, and K contents are commonly measured. For example, the GSNL till geochemistry dataset yielded over 17 000 points while the National Soil Network yielded 231 samples in Newfoundland, only 65 of which recorded Al values. Further, when pedogenesis originates from glacial till, soil forms a direct and often indistinguishable boundary with the till. Till samples thus provide an excellent alternative to bedrock maps for soil parent material.

4.2 Novelty of this work

While previous mapping efforts have addressed the distribution of reactive Al in soils (Ren et al., 2024), this is the first attempt to map it in soil parent material. Support for modelling till geochemistry through existing datasets is likely more common than reactive Al in soil, as demonstrated in Newfoundland, with potential for providing a key factor supporting SOC modelling across more regions. Few studies investigate till geochemistry from a predictive geochemical standpoint, and those that do focus on economically viable trace elements such as rare earth elements (e.g. Casey et al., 2024). The results developed from this model can be used directly in future SOC models in Newfoundland and provides an important demonstration relevant to the development of SOC modelling in other boreal regions.

The uncertainty calibration is an uncommon yet important component to this project. Addressing uncertainty in SOC modelling has been identified as essential (Mishra et al., 2013); however, while many digital soil mapping studies produce spatial predictions (Pouladi et al., 2023) few provide spatial estimates of uncertainty, and even fewer provide calibrated uncertainty estimates. The calibration step may be critical for accurately reporting model uncertainty; here, calibration improved uncertainty estimates by over 30 % in both cases. These calibrated estimates of uncertainty are critical for operationalizing these results to accurately constrain estimates of SOC stocks. An accurate spatial representation of uncertainty indicates regions of the study area that are predicted with lower confidence, which could be (a) integrated into subsequent models or calculations of SOC stock for the island of Newfoundland, or (b) used to guide future sampling to improve spatial estimates of Al %.

4.3 Limitations

Several constraints limit these models. While our results indicate bedrock geology is the largest control on Al % generally, there is a lack of understanding as to what variables, or interactions between variables, exactly control Al weathering within glacial till. For example, biotic factors, which have been shown to influence soil Al % at coarser scales (Ren et al., 2024), were not explored for the purpose of this study. Further, the model used geological context to leverage spatial autocorrelation within the samples, mimicking a geostatistical approach. While this improved model performance, it renders partial dependence plots for the spatial model difficult to interpret, and the individual variables difficult to rank. Finally, because the work was focused strongly on spatial predictive accuracy, multicollinearity was left largely unaddressed, making it possible that variable importance was split across multiple correlated predictor variables. Such limitations mean that this modelling approach is less useful for discerning strong predictors for available Al content in glacial till but rather is a useful approach for obtaining a spatial prediction for available Al in soil parent material. This represents a first order component for predictive modelling of boreal forest SOC. Next steps, however, should focus on assessing predictors of the relationship between available Al in parent material and reactive Al in the overlying boreal forest mineral soil including biotic topographic, and climatic factors and their interactions (Brantley et al., 2017; Maher, 2010, Patrick 2023). Finally, the CIA is a one-dimensional value quantifying the proportion of mobile cations to immobile Al, and is the most common index relying on field and geochemical studies, verified by kinetic and thermodynamic modeling decades ago (Nesbitt and Young, 1982, 1984; Nesbitt et al., 1997). Another index using commonly available major element geochemistry, such as the mafic index of alteration (MIA; Babechuk et al., 2014), could also be applied. However, all such one-dimensional indices suffer from mineralogical controls independent of petrogenic effects (Nesbitt and Wilson, 1992; Babechuk et al., 2014). These are issues largely solved for the purposes of siliciclastic sedimentary geochemistry using three-dimensional approaches (Babechuk and Fedo, 2023). Furthermore, use of the chemical index of alteration (CIA) to inform on available Al in till is likely to vary with large variation in carbonate content. If incongruencies between in situ soil reactive Al and Al % are found, then use of complementary data such as total inorganic carbon could be evaluated as a means to improve the relationship between Al availability in the till and overlying mineral soil.

4.4 Utility of the spatial and aspatial models

The patterns of uncertainty across the spatial and aspatial models indicate that the spatial model may be more useful in previously sampled areas, while the aspatial model may deliver higher certainty in areas with scarce samples. While the overall accuracy of the aspatial model is lower, uncertainty is distributed more evenly across space than the spatial model. As seen in Figs. 4 and 5 the calibration accuracy is slightly higher in the spatial model than the aspatial model. It may be possible to aggregate both models for regional scale mapping, while for finer scale mapping the model with higher local performance will be preferred.

4.5 Future work

Soil geochemistry is controlled largely by parent material – in this case till – and this in turn controls the ability of soils formed to store C. For example, reactive Al content of soil is a key factor controlling SOC content of boreal forest podzols in the study region (Patrick et al., 2022). The exact nature of the relationship between till and soil geochemistry is not fully understood, however, topography and climate likely play key interacting roles. To inform this relationship, soil and till should be sampled simultaneously at sites informed by this predictive model. Specifically, sampling where there is forecasted high or very low Al, as these extremes may best illustrate the relationships between soil Al % and till Al %. Sampling should be conducted on similar slopes to ensure that the relationship between till and soil is being captured without influence from terrain characteristics, where higher slope areas experience less vertical and more lateral flow resulting in a shallower horizon of SOC enrichment (Patrick, 2023). As depth of horizon is strongly correlated to other soil characteristics governing SOC storage and this may cause variations in SOC content not attributable to parent material alone (Patrick, 2023). Sampling on high and low slopes should then also proceed to assess controls of terrain on SOC storage.

The GSNL continually updates their till sample geochemistry database and new samples should be incorporated as they become available. Till samples from western Newfoundland, particularly the Great Northern Peninsula, will be especially valuable. Additionally, independent glacial till sampling should be prioritized in high uncertainty areas to increase the robustness of these models.

As outlined in Sect. 2.3, Al % is calculated using whole rock elemental composition, and therefore is a proxy for the source of Al to overlying soils rather than a direct estimate of the reactive Al pool contributing to OMC formation in soil. This reflects the aim of this study: to develop an approach for modelling a proxy for a key geological factor influencing SOC sequestration. Given the prevalence of till geochemistry datasets, this is a feasible approach compared to the use of reactive metal content in soils. This does mean it must be tested as an input to soil C models though. Future work should include selective dissolution extraction to directly measure reactive Al in soils, followed by a comparison of this to underlying till Al % values to elucidate the relationship between proxy and soil pool.

Given the important role that parent material plays in regulating SOC (Slessarev et al., 2022; Patrick et al., 2022) increased attention must be given to improving parent material composition estimates for SOC modelling. This includes further investigation into key parameters controlling Al weathering and thus reactive Al in parent material. For example, integrating remotely sensed terrain covariates such as slope aspect and topographic wetness index into future modelling of Al % may strengthen predictions.

The results of the calibration in this study demonstrate the advantages of uncertainty quantification. Future SOC modelling should incorporate this dimension, as should any future mapping of parent material, as it is crucial to producing meaningful model predictions. In the case of modelling SOC, the more reliable the uncertainty the clearer the utility of the predictions for future model development or application for mitigation efforts.

5 Conclusion

Using both a spatial and aspatial approach, available Al % was successfully modelled and mapped in soil parent material across the Island of Newfoundland for the first time at a high resolution. Calibrated uncertainty was also estimated, providing spatial estimates of model confidence throughout the study area. There is broad utility in predicting parent material geochemistry because it is known to influence soil properties and modelling glacial till draws on a sample base much wider, in many cases, than available soil samples. While future research may expand the sample distribution, and explore relationship between glacial till and soil geochemistry, this work provides a critical component for a future mechanistic SOC model tailored to the soils of Newfoundland and Labrador. This will enable further investigation of the performance of previous modelling across various landscapes, topographies, and soil Al % delivering valuable insight into the spatial distribution of boreal forest SOC across the island. By enabling continuous predictions of SOC and better estimates of the terrestrial carbon reservoir in Newfoundland, this future work will provide a template for such research in other boreal regions.