Investigating the Spatial Variability in Soil Geochemical and Colour Properties Across Two Contrasting Land Uses in South-Central Manitoba

1 Introduction

Variation in soil biological, chemical, and physical properties occurs across the landscape and in response to both regional and local (i.e., field-scale) variations in the five soil forming factors: parent material, relief or topography, biota, climate, and time. Superimposed on this is the influence of changes in land use and current and historic management practices which can further modify soil properties. Quantifying and understanding the patterns and drivers for this variation is an important component of many agri-environmental studies. For example, to meet the desired level of precision for agronomic and environmental nutrient management plans the spatial variability in soil nutrients will influence the soil sampling design in terms of number and locations of soil samples (Kariuki et al., 2009; Starr et al., 1995).

Sediment source fingerprinting is a watershed-scale technique that is used to identify and quantify the relative proportions of sediment derived from unique sources. This technique uses natural occurring biogeochemical properties as fingerprints (i.e., tracers) to discriminate between potential sources of sediment and are linked to downstream sediment using mixing models. From a sediment fingerprinting perspective, investigating the spatial variability of soil properties at a watershed-scale can be advantageous to identify, classify, and distinguish between potential sources of sediment (Pulley et al., 2017). However, investigating spatial variability at smaller scales is less common (Collins et al., 2019; e.g., Du and Walling, 2017; Luna Miño et al., 2024; Pulley et al., 2018) and remains a research priority (Collins et al., 2020).

There are three main, interconnected, ways that spatial variability in fingerprint properties are an important aspect of sediment fingerprinting. First is to adequately quantify the fingerprint properties such that it is representative of that source. For some fingerprints the variability is not random but rather varies in a more systematic way. For example, the pattern of fallout radionuclides will reflect the patterns of soil erosion and deposition (Wilkinson et al., 2015). Designing and implementing source sampling plans need to take this into consideration as the sampling designed used has been shown to influence the characterization of wide range of commonly used fingerprints (Luna Miño et al., 2024).

Secondly, the issue of spatial variability of fingerprint properties is further complicated by overlying spatial variability in the rates of erosion and sediment delivery. Incorporation of both types of variability into the mixing model will provide a more reliable estimate of the proportion of sediment derived from each source. Many mixing models have well defined inputs (sources) and outputs (sediment) that are characterized by their mean and standard deviation and the spatial distribution or pattern of fingerprints are not considered. This is not ideal as the values of samples that are collected closer, and more hydrologically connected, to the stream network may in fact present a better representation of that source despite potentially deviating from the mean value. This issue can be addressed by strategic sampling where the more likely to erode areas are targeted for sampling. However, a considerable amount of information and insight is lost through that approach. There has been some progress using information on erosion rates to calculate a erosion rate-weighted mean (Du and Walling, 2017; Wilkinson et al., 2015) and using spatially interpolated maps of fingerprint values to provide a finer resolution of the fingerprint variability within each source (Haddadchi et al., 2019).

Lastly, understanding the geomorphic, hydrologic, and biochemical processes that have led to the observed patterns in spatial variability helps in the selection of robust and reliable fingerprints and/or guide the sampling design for source characterization. In selecting fingerprints that provide good discrimination between sources many studies typically used a statistical-based approach (Collins et al., 1997). However, there are concerns that this approach may result in the inclusion of false positives (i.e., type I error) or non-conservative fingerprints (Koiter et al., 2013a). Consequently, there has been a call for the inclusion of a process-based (e.g., weathering, erosion) or geologic/lithologic-based explanation of the fingerprints selected to address these concerns (Collins et al., 2020). Furthermore, there is also a lack of standardization in how sediment source areas are sampled (e.g., judgement, random, transect, grid, stratified) and it can be difficult to have an efficient sampling design without prior knowledge of why and how soil properties vary across the landscape (Luna Miño et al., 2024). Prior knowledge of the spatial variability of soil fingerprint properties would be beneficial; however, in practice this can be difficult, particularly with geochemical properties as routine lab analysis often return information on more than 50 elements. The spatial patterns of some soil properties are well studied because of their agronomic importance or ability to infer other important soil properties and processes and can include fallout radionuclides [e.g., ¹³⁷Cs, ⁷Be; (Ritchie et al., 1970)], plant nutrients [e.g., N, P; (Vasu et al., 2017)], soil colour [e.g., hue, value; (Viscarra Rossel et al., 2006)], major non-acid forming cations [e.g., Ca, Na; (Sun et al., 2021a)]. In contrast, the processes that control the distribution of other soil properties, such as rare earth elements and trace metals, are less well studied or tend to be site-specific, making it difficult to draw generalizations.

Terrain attributes such as elevation, slope curvature, slope position, and soil wetness indices have been shown to be useful information in the understanding and modelling of a range of soil properties including soil moisture (Beaudette et al., 2013), texture (Kokulan et al., 2018), colour (Brown et al., 2004), organic matter (Zhang et al., 2012), conductivity (Umali et al., 2012), and geochemistry (Lima et al., 2023). Similar techniques may provide additional insight into the pedologic and geomorphic processes that drive the observed patterns of fingerprint properties within a given source. Since digital elevation models (DEMs) are more publicly available and can in some case be generated using drone imagery, while soil property data are often limited, terrain attributes derived from DEMs can be used to guide sampling design.

This study builds on the previous work of Luna Miño (Luna Miño et al., 2024) where the impact of three different sampling designs on the characterization of source materials, within the framework of the sediment fingerprinting approach, was assessed. This study expands that study by using the data from grid sampling approach to assess the spatial autocorrelation, create soil property (i.e., fingerprint) maps, and identify important terrain attributes driving the observed patterns. The objectives of this study were (1) to investigate the spatial variability of a range of soil colour and geochemical properties in an agricultural and forested site; and (2) to assess the relative importance and correlation of terrain attributes with the spatial distribution of these soil properties. Together, these objectives address how terrain attributes may be used to understand spatial distributions of soil properties and help guide sampling design.

2 Methods

2.1 Site description

Two sites of contrasting land uses located in the Wilson Creek Watershed (WCW), near McCreary, Manitoba, Canada were selected to investigate the spatial variability in fingerprint properties. The headwaters of the WCW are located on top of the Manitoba Escarpment within the boundary of Riding Mountain National Park. There is a ~300m drop in elevation crosses the escarpment where the streams become deeply incised. At the base of the escarpment is a large alluvial fan situated in the lacustrine deposits of glacial lake Aggasiz where the main stem has a meandering form. However, beyond the national park boundary the stream flows straight through an engineered drain until it reaches the Turtle river (Figure 1). Both sites are both hydrologicaly connected to the mainstem of the Wilson Creek

The first site was a mixedwood forest including white and black spruce (Picea glauca, Picea mariana), balsam fir (Abies balsamea), larch (Larix laricina) and young stands of deciduous trees including trembling aspen (Populus tremuloides). The forested site is located within the boundaries of the national park where there is little disturbance beyond recreational hiking trails. The soils within the park are not well mapped but likely are part of the Grey Wooded soil association (Luvisol) consisting of fine sandy loam to clay loam soils developed on boulder till of mostly shale with some limestone, and granitic rocks (Ehrlich et al., 1958). The second site is under agricultural production and includes rotations of grain crops and forage. The site is mapped to the Edwards Soil Series (Cumulic Regosol) consisting of silty clay loam to silty clay soil developed on recent alluvial deposits (Ehrlich et al., 1958).

The Köppen-Geiger climate classification of the WCW is cold, without dry season, and with warm summer (Dfb) (Beck et al., 2018). The average annual precipitation is ~539 mm, with approximately 27% falling as snow with a mean annual temperature is 3.0°C (Environment and Canada, 2024). The hydrology of the watershed is snowmelt dominated with ~ 80% of the cumulative runoff occurring during the spring season (May and June) (MacKay, 1970).

In [1]:

p6 <- p1 + guide_area() + p2 + p3 + plot_layout(ncol = 2, guides = "collect") + plot_annotation(tag_levels = 'a', tag_suffix = ')')
p6

Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
give correct results for longitude/latitude data
Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
give correct results for longitude/latitude data
Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
give correct results for longitude/latitude data
Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
give correct results for longitude/latitude data
Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
give correct results for longitude/latitude data

Figure 1: a) Map showing the location of the study sites within Canada. Location of the two study sites and nearby town of McCreary, and regional b) land use, and c) topography.

2.2 Soil sampling and analysis

This study uses samples and data collected as part of the grid sampling design outlined in Luna Miño (Luna Miño et al., 2024). Briefly, at each site 49 samples were collected using a soil auger on a 7x7 grid at a 100m spacing. Within the forested surface soil samples were collected below the LFH layer to a depth of 5cm, and the agricultural site was sampled to a depth of 15cm to account for the regular mixing of the soil due to tillage and other field operations.

Samples were dried, homogenized with a mortar and pestle, and sieved through a 63 𝜇m sieve to remove the sand fraction. The sand fraction was removed in an effort to reduce the differences in grain size and organic matter content between the two sites (Laceby et al., 2017). Samples were analyzed for a broad suite geochemical element using inductively coupled plasma mass spectrometry (ICP-MS) following a microwave-assisted digestion with aqua-regia (ALS Mineral Division, North Vancouver, BC, Canada). Spectral measurements were made with a spectroradiometer (ASD FieldSpecPro Malvern Panalytical Inc Westborough MA 01581, United States). Spectral reflectance measurements were taken in 1 nm increments over the 0.4-2.5 μm wavelength range. Both samples and Spectralon standard (white reference) were illuminated with a white light source using a halogen-based lamp (12 VDC, 20 Watt). Following the method outlined in Boudreault et al. (Boudreault et al., 2018), fifteen colour coefficients (R, G, B, x, y, Y, X, Z, L, a*, b*, u*, v*, c*, h*) were calculated for each sample (Koiter, 2021). Based on the results of Luna Miño (Luna Miño et al., 2024), a composite fingerprint consisting of 10 geochemical elements (Ca, Co, Cs, Fe, Li, La, Nb, Ni, Rb, and Sr) and five colour coefficients (a*, b*, c*, h*, and x) were identifying as providing a strong discrimination between the agricultural and forested surface soils. These fifteen soil properties are the focus of the detailed spatial analysis detailed in this study.

2.3 Geospatial and terrain analysis

All geostatistics were performed with ArcGIS Pro (v 3.3.0 Esri, 2024). Semivariograms were used to quantify spatial correlation for each of the 15 soil properties. The optimization tool, based on minimizing the mean square error, was used to parameterize the semivariogram model. Kriging was used to interpolate and generate maps of each soil property. The exploratory interpolation tool (Geostatistical Analyst extension) was used to select the kriging type with the highest ranked prediction accuracy.

A Digital Elevation Model (DEM) for the forested site was acquired from publicly available data (Canada, 2024). A DEM for the agricultural site was generated by photogrammetry using UAV imagery, including the use of ground control and check points, with Agisoft Metashape Professional (v1.8.2 Agisoft, 2021). Ordinary kriging was used to calculate a 1 m gridded digital elevation model for each site. Geographic information software (SAGA v2.1.4 Conrad et al., 2015) was used to calculate six additional terrain attributes and included plan and profile curvatures, saga wetness index, catchment area, relative slope position, and vertical channel network distance (Table S2).

2.4 Data analysis

All subsequent statistical analysis was undertaken using R statistical Software v4.4.0 (Team, 2024) through RStudio Integrated Development Environment v2024.04.2 (RStudio, 2024). Plots and maps were created using the R package ggplot2 v 3.5.2 (Wickham, 2016). Skewness was categorized as values between -0.5 and 0.5 considered approximately symmetric, -1.0 to -0.5 or 0.5 to 1 as moderately skewed, and < -1.0 or > 1.0 as highly skewed. Coefficient of variation (CV) thresholds were categorized as low (<15%), moderate (15–35%), high (35–75%), and very high (>75%). Interpolated soil property and terrain attribute data were resampled to a 10 m resolution prior to analysis (terra v1.8.29 Hijmans, 2024). Random Forest Regression (randomForest v4.7.1.2 Liaw and Wiener, 2002) was used to assess the relative importance of the terrain attributes on the spatial distribution of soil properties. The dataset was randomly split into training, validation, and testing datasets. Multicollinearity among the terrain attributed was assessed using the Variance Inflation Factor with a threshold of eight and correlated terrain attributes were removed (usdm v2.1.7 Naimi et al., 2014). The number of variables randomly sampled as candidates at each split within the random forest model was tuned using the training and validation data sets (caret v7.0.1 Kuhn and Max, 2008). The number of trees to grow was set to 500 and model performance was assessed using the Mean Square Error (MSE) and percent variance explained for both the training (Out of Bag Error) and the validation data sets. To test the model, actual and predicted values were plotted and the R^₂ and MSE were calculated using the testing data set. Because analyzing interpolated data can cause issues, the random forest model was used to predict the original 49 non-interpolated data points at each site as an additional check.

3 Results

3.1 Univariate summary

Overall, the agricultural site had soil colour and geochemical properties that exhibited lower variability and more symmetrical data distributions as compared to the forested site (Table 1). All 15 colour properties at both sites displayed approximately symmetrical distributions. At the agricultural site, all colour properties were characterized by low coefficients of variation (CV), while the forested site showed slightly greater variability, with 10 colour properties having low CVs and five having moderate CVs.

Similarly, geochemical data at the agricultural site showed lower variability and greater symmetry. Most elements were approximately symmetrical, with only nine exhibiting moderate skewness and five highly skewed (Table 1). Variability was also limited, with the majority of elements having low CVs; 12 had moderate CVs and five had high CVs. In contrast, the forested site showed greater skewness and variability: seven elements exhibited moderate skewness, 14 were highly skewed, 28 had moderate CVs, six had high CVs, and two had very high CVs.

In [2]:

In [3]:

summary_tab |>
  group_by(site) |>
  select(-group, -se) |>
  gt() |>
  #fmt_markdown(columns = Element)|>
  fmt_number(decimals = 2) |>
  #cols_align(align = "left") |>
  tab_style(style =  cell_text(weight = "bold", align = "center"), locations =  cells_row_groups()) |>
  tab_options(column_labels.font.weight = "bold") |>
tab_style(
    style = list(cell_text(style = "italic")),
    locations = cells_body(
      columns = Property,
      rows = Property %in% c("a*", "b*", "c*", "h*", "x")))

Table 1: Summary univariate statistics of selected geochemical and colour soil properties for each site (n = 49). Geochemical concentrations are reported in ppm, excecpt Ca and Fe(%).

Property	Mean	SD	Max	Min	Skewness	CV
Agriculture
Ca	4.00	2.19	8.78	0.95	0.28	54.66
Co	8.76	0.83	10.60	7.50	0.52	9.48
Cs	0.75	0.15	1.07	0.47	0.18	19.93
Fe	1.92	0.09	2.11	1.71	−0.25	4.70
Li	15.62	1.42	19.80	12.80	0.62	9.11
La	18.23	1.22	20.20	15.50	−0.29	6.71
Nb	0.59	0.06	0.73	0.46	0.45	9.67
Ni	29.63	2.72	35.70	25.00	0.36	9.17
Rb	18.43	4.33	26.70	10.20	0.24	23.48
Sr	91.31	38.98	163.50	38.60	0.09	42.69
a*	3.38	0.32	4.15	2.59	−0.03	9.53
b*	8.84	0.97	10.59	6.69	−0.18	11.00
c*	9.47	1.02	11.32	7.17	−0.19	10.74
h*	1.20	0.01	1.23	1.18	0.19	1.12
x	0.47	0.00	0.48	0.47	0.06	0.46
Forest
Ca	1.89	1.53	5.46	0.47	1.07	81.12
Co	6.76	1.39	9.60	4.00	0.03	20.62
Cs	0.55	0.12	0.78	0.34	0.25	21.73
Fe	1.18	0.13	1.46	0.83	−0.58	11.24
Li	6.47	0.90	8.60	4.30	−0.02	13.89
La	15.00	2.60	21.80	10.30	0.33	17.31
Nb	0.37	0.06	0.56	0.17	−0.68	17.10
Ni	18.09	3.90	28.00	11.00	0.33	21.55
Rb	13.83	1.85	18.10	9.90	0.27	13.40
Sr	32.43	12.60	64.20	15.30	0.98	38.87
a*	5.73	0.41	6.56	4.41	−0.38	7.10
b*	12.47	2.01	15.91	8.02	0.22	16.11
c*	13.74	1.94	17.00	9.15	0.15	14.15
h*	1.13	0.05	1.23	1.06	0.34	4.13
x	0.49	0.00	0.49	0.48	−0.21	0.47

The agricultural site has a relatively flat topography with an elevation change of approximately 3m, with the field draining toward a ditch in the northeast corner. The forested site has a relatively more complex topography, with a channel flowing from the southwest toward the northeast and an overall elevation difference of 18 m across the site. The mean plan and profile curvature measurements for both sites are near zero indicating an area of sediment transit and not accumulation or erosion (Table 2). The agricultural site had a higher SAGA Wetness Index but the forested site had a larger range in values and exhibited a higher degree of variability. The forested site exhibited a smaller mean Relative Slope Position value (streams and depressional areas) and a smaller Vertical Distance to Channel Network, and for both terrain attributes a greater variability as compared to the agricultural reflecting the presence of the stream crossing the forested site.

In [4]:

In [5]:

terrain_summary |>
  rename("Property" = "property") |>
  group_by(site) |>
  gt() |>
  #fmt_markdown(columns = Property)|>
  fmt_number(n_sigfig = 3) |>
  fmt_scientific(rows = Property %in% c("Plan Curvature", "Profile Curvature", "Vert. Dist. Channel")) |>
  #fmt_number(row = everything() > 100, decimal = 0) |>
  #cols_align(align = "left") |>
  tab_style(style =  cell_text(weight = "bold", align = "center"), locations =  cells_row_groups()) |>
  tab_options(column_labels.font.weight = "bold") |>
tab_style(
    style = list(cell_text(style = "italic")),
    locations = cells_body(
      columns = Property,
      rows = Property %in% c("a*", "b*", "c*", "h*", "x")))

Table 2: Summary statistics for the interpolated values (10m resolution) for selected geochemical and colour soil properties and terrain attributes for each site. Geochemical concentrations are reported in ppm, excecpt Ca and Fe(%).

Property	Mean	SD	Max	Min	Skewness	CV
Agriculture
Ca	4.12	2.10	8.76	0.918	0.0727	51.0
Co	8.75	0.664	10.6	7.52	0.431	7.59
Cs	0.729	0.123	1.07	0.458	0.376	16.9
Fe	1.92	0.0644	2.10	1.73	−0.450	3.36
Li	15.7	1.16	19.3	13.2	0.551	7.38
La	18.2	0.817	19.8	16.5	−0.268	4.49
Nb	0.593	0.0550	0.740	0.459	0.569	9.27
Ni	29.9	2.23	34.5	26.3	−0.0100	7.46
Rb	18.0	3.94	26.1	11.5	0.498	21.8
Sr	93.4	38.6	167	36.3	0.00105	41.3
a*	3.34	0.211	3.83	2.88	0.0621	6.33
b*	8.73	0.707	10.2	6.98	−0.162	8.10
c*	9.34	0.762	11.0	7.41	−0.158	8.15
h*	1.20	0.00977	1.23	1.18	−0.0603	0.811
x	0.473	0.00149	0.477	0.470	−0.0168	0.314
Plan Curvature	1.65 × 10−6	1.36 × 10−4	6.57 × 10−4	−5.07 × 10−4	3.54 × 10−1	8.24 × 103
Profile Curvature	−7.64 × 10−6	1.53 × 10−4	5.83 × 10−4	−6.47 × 10−4	9.51 × 10−2	−2.00 × 103
SAGA Wetness Index	9.64	0.704	11.2	7.77	−0.122	7.30
Catchment Area	475	1,010	10,100	4.35	4.76	213
Rel. Slope Position	0.718	0.288	1.20	0.0221	−0.946	40.1
Vert. Dist. Channel	5.98 × 10−2	4.10 × 10−2	2.92 × 10−1	4.25 × 10−3	1.21	6.85 × 101
Elevation	310	0.593	312	309	0.615	0.191
Forest
Ca	1.88	0.769	3.61	0.787	0.202	40.8
Co	6.80	0.632	8.66	4.93	−0.200	9.30
Cs	0.551	0.0737	0.714	0.423	0.297	13.4
Li	6.43	0.694	8.46	4.39	−0.136	10.8
La	15.0	1.57	18.5	11.5	−0.0324	10.4
Nb	0.370	0.0356	0.440	0.278	−0.436	9.64
Ni	18.2	2.49	24.9	14.3	0.314	13.7
Sr	31.6	8.50	53.1	18.1	0.716	26.9
h*	1.13	0.0371	1.22	1.06	0.257	3.27
Plan Curvature	3.97 × 10−4	3.27 × 10−3	2.89 × 10−2	−2.62 × 10−2	7.91 × 10−1	8.22 × 102
Profile Curvature	−1.83 × 10−4	9.47 × 10−3	6.37 × 10−2	−7.37 × 10−2	−5.31 × 10−1	−5.18 × 103
SAGA Wetness Index	6.00	0.988	8.48	2.21	−0.430	16.5
Catchment Area	571	1,940	25,400	3.44	6.60	339
Rel. Slope Position	0.222	0.232	0.993	0.00617	1.56	105
Vert. Dist. Channel	4.15 × 10−1	4.43 × 10−1	3.66	2.02 × 10−2	2.96	1.07 × 102
Elevation	369	3.34	377	359	−0.184	0.904

3.2 Spatial analysis

Soil colour and geochemical composition varied across both sites. In the agricultural field, all 15 soil color and geochemical properties exhibited spatial autocorrelation with most properties demonstrating a strong spatial dependency (Table 3). Some of the soil properties presented a pattern that roughly matches (e.g., Rb, Cs) or mirrors (e.g., Ca, Sr) the overall topography of the site with a gradation between the highest point in the south-west corner towards the lowest points in the north-east (Figure 2). Other properties appear to have more localized highs and low concentrations/values (e.g., c*, h). The geochemical concentrations of Ca and Rb had the largest range values and as result displayed a less patchy distribution across the site. The nugget (Co) was small for all soil properties (<1.5), and Sr had an exceptionally large sill value (900).

At the forested site, the geochemical concentrations of Fe and Rb, along with the color properties a*, b*, c*, and x showed no spatial autocorrelation and were excluded from further analysis and four and five properties exhibiting strong and moderate spatial dependency, respectively (Table 3). In comparison to the agricultural site, the soil properties at the forested site displayed a more moderate spatial dependency. The nugget (Co) was generally small for most soil properties (<2) with the exception of La and Ni. The range values were similar across the different soil properties and fell between 176 and 298 m. Overall, the influence of the channel and floodplain environment can be seen in the pattern of the nine soil properties (Figure 3).

In [6]:

In [7]:

vario |>
  group_by(site) |>
  gt() |>
  fmt_number(decimal = 1) |>
  fmt_number(columns = c("Range (m)", "C/(C + Co) (%)"), decimal = 0) |>
  cols_label(r2 = "{{r^2}}") |>
  tab_footnote(
    footnote = "Models are all isotropic.",
    locations = cells_column_labels(columns = "Kriging Type"))  |>
  tab_footnote(
    footnote = "Strong spatial dependency (C/(C + Co) % >75); Moderate spatial dependency (C/(C + Co) % between 75 and 25); Low spatial dependency (C/(C + Co) % <25).",
    locations = cells_column_labels(columns = `Spatial Class`))  |>
  tab_style(style =  cell_text(weight = "bold", align = "center"), locations =  cells_row_groups()) |>
  tab_options(column_labels.font.weight = "bold")

Table 3: Geostatistical parameters of the fitted semivariogram models of selected colour and geochemical properties within the agricultural and forested sites.

Property	Kriging Type1	Nugget (Co)	Sill (Co + C)	C/(C + Co) (%)	Range (m)	r2	Spatial Class2
Agriculture
Ca	Universal	0.0	7.2	100	580	0.9	Strong
Co	Simple	0.0	0.7	100	208	0.4	Strong
Cs	Ordinary	0.0	0.0	100	210	0.5	Strong
Fe	Ordinary	0.0	0.0	100	185	0.2	Strong
Li	Universal	0.3	1.5	81	185	0.6	Strong
La	Simple	0.4	1.0	56	308	0.5	Moderate
Nb	Universal	0.0	0.0	91	210	0.7	Strong
Ni	Ordinary	1.4	8.9	84	352	0.6	Strong
Rb	Ordinary	1.4	27.6	95	551	0.9	Strong
Sr	Ordinary	0.9	900.2	100	220	1.0	Strong
a*	Ordinary	0.4	1.0	59	288	0.3	Moderate
b*	Simple	0.2	0.9	83	199	0.3	Strong
c*	Simple	0.1	0.9	87	199	0.3	Strong
h*	Simple	0.0	1.1	100	185	0.2	Strong
x	Simple	0.4	1.0	58	220	0.1	Moderate
Forest
Ca	Ordinary	1.6	2.7	41	269	0.2	Moderate
Co	Ordinary	0.0	2.1	100	298	0.1	Strong
Cs	Ordinary	0.0	0.0	83	237	0.2	Strong
Li	Ordinary	0.0	0.8	100	222	0.3	Strong
La	Ordinary	3.1	7.4	59	176	0.1	Moderate
Nb	Ordinary	0.0	0.0	51	224	0.2	Moderate
Ni	Universal	6.7	15.8	57	187	0.2	Moderate
Sr	Simple	0.4	1.0	65	229	0.4	Moderate
h*	Universal	0.0	0.0	100	230	0.3	Strong
1 Models are all isotropic.
2 Strong spatial dependency (C/(C + Co) % >75); Moderate spatial dependency (C/(C + Co) % between 75 and 25); Low spatial dependency (C/(C + Co) % <25).

In [8]:

p_ag <- ca + co + cs + fe + li + la + nb + ni + rb + sr + a_col + b_col + c_col + h_col + x_col + dem + plot_layout(ncol = 4) 

p_ag

#ggsave(filename = "../Submission/Figure 2.png", plot = p_ag,  height = 240, width = 190, units = "mm", dpi = 600)

Figure 2: Kriged maps of select colour and geochemical properties and elevation across the agricultural site.

In [9]:

p_for <- ca_forest + co_forest + cs_forest + la_forest + li_forest + nb_forest + ni_forest + sr_forest + h_col_forest + dem_forest +
  plot_layout(ncol = 4) 
p_for

#ggsave(filename = "../Submission/Figure 3.png", plot = p_for,  height = 175, width = 190, units = "mm", dpi = 600)

Figure 3: Kriged map of select colour and geochemical properties and elevation across the forested site.

Across both sites, there was a significant (p < 0.05) correlation between the selected soil properties and the terrain attributes, with the exception of the plan and profile curvature attributes (Table S3). The elevation attribute generally had higher correlation coefficients; however, the direction and strength of the correlation did vary between both site and soil property. Overall, the random forest regression models exhibited relatively strong predictive performance, with the models better performing at the agricultural site compared to the forested site (Table 4). With the exception of the Ni concentration and x colour values at the agricultural site, elevation was consistently the terrain attribute that provided the greatest predictive power (Figure 4). SAGA Wetness and relative slope position were generally the second and third most informative terrain attributes. Plan curvature was consistently ranked least important predictive terrain attribute.

In [10]:

In [11]:

m_performance |>
  group_by(site) |>
  gt() |>
  tab_spanner(label = "Training",
               columns = c(Var_exp, MSE)) |>
  tab_spanner(label = "Validation",
              columns = c(Var_exp_test, MSE_test)) |>
  tab_spanner(label = "Test",
              columns = c(R2, mse_test)) |>
    tab_spanner(label = "Original",
              columns = c(R2_org, MSE_org)) |>
  cols_move_to_start(columns = property) |>
  # tab_spanner(label = "Training", columns = c("MSE", "Var_exp")) |>
  # tab_spanner(label = "Validation", columns = c("MSE_test", "Var_exp_test")) |>
  # tab_spanner(label = "Testing", columns = "R2") |>
  fmt_number(columns = c("MSE_test", "MSE", "mse_test", "MSE_org"), decimal = 2) |>
  fmt_number(columns = c("Var_exp_test", "Var_exp"), decimal = 1) |>
  fmt_number(columns = c("R2", "R2_org"), decimal = 2) |>
  tab_style(style =  cell_text(weight = "bold", align = "center"), locations =  cells_row_groups()) |>
  tab_style(style =  cell_text(align = "left"), locations =  cells_body()) |>
  tab_style(style =  cell_text(weight = "bold", align = "left"), locations =  cells_column_labels()) |>
  # cols_label(Var_exp = md("% Var<br/>Training"), Var_exp_test = md("% Var<br/>Validation"), MSE_test = md("MSE<br/>Testing"), MSE = md("MSE<br/>Training"), R2 = md("{{R^2}}<br/>Testing"), property = "Property") |>
    cols_label(Var_exp = "% Variance", Var_exp_test = "% Variance", MSE_test = "MSE", MSE = "MSE", R2 = "{{R^2}}", R2_org = "{{R^2}}", property = "Property", mse_test = "MSE", MSE_org = "MSE") |>
  row_group_order(groups = c("Agriculture", "Forest"))  |>
  tab_style(
    style = list(cell_text(style = "italic")),
    locations = cells_body(
      columns = property,
      rows = property %in% c("a*", "b*", "c*", "h*", "x"))) |>
    tab_footnote(
    footnote = "Mean square error",
    locations = cells_column_labels(columns = c("mse_test", "MSE_test", "MSE", "MSE_org"))) |>
  tab_footnote(
    footnote = "Percent variance explained",
    locations = cells_column_labels(columns = c("Var_exp", "Var_exp_test"))) 

Table 4: Model summary and performance statistics for the random forest regression using the training, validation, test and original (non-interpolated) data sets.

Property	Training		Validation		Test		Original
	% Variance1	MSE2	% Variance1	MSE2	R2	MSE2	R2	MSE2
Agriculture
Ca	91.6	0.37	91.8	0.36	0.91	0.38	0.95	0.23
Co	79.8	0.09	82.5	0.08	0.80	0.08	0.88	0.08
Cs	85.7	0.00	86.4	0.00	0.85	0.00	0.92	0.00
Fe	69.6	0.00	70.9	0.00	0.69	0.00	0.83	0.00
Li	59.3	0.54	59.8	0.53	0.64	0.51	0.88	0.24
La	93.0	0.05	93.1	0.04	0.93	0.05	0.96	0.03
Nb	57.3	0.00	59.1	0.00	0.55	0.00	0.71	0.00
Ni	93.1	0.34	93.7	0.33	0.93	0.34	0.95	0.25
Rb	95.3	0.73	96.1	0.64	0.95	0.79	0.98	0.39
Sr	93.5	97.22	93.6	93.97	0.93	105.59	0.97	44.77
a*	85.0	0.01	86.9	0.01	0.85	0.01	0.91	0.00
b*	72.5	0.14	75.3	0.12	0.72	0.15	0.89	0.09
c*	73.2	0.15	75.9	0.14	0.73	0.17	0.89	0.10
h*	58.3	0.00	58.6	0.00	0.56	0.00	0.73	0.00
x	73.3	0.00	73.6	0.00	0.69	0.00	0.82	0.00
Forest
Co	39.1	0.24	42.9	0.23	0.48	0.21	0.77	0.29
Cs	64.1	0.00	67.1	0.00	0.66	0.00	0.86	0.00
Li	41.3	0.28	42.0	0.28	0.46	0.28	0.66	0.26
La	43.3	1.40	47.5	1.32	0.48	1.23	0.78	0.60
Nb	55.0	0.00	55.9	0.00	0.58	0.00	0.84	0.00
Sr	59.4	29.43	59.1	29.66	0.59	29.25	0.82	13.88
h*	58.8	0.00	60.3	0.00	0.62	0.00	0.86	0.00
1 Percent variance explained
2 Mean square error

In [12]:

p <- g2 + g1 + plot_layout(ncol = 1, heights = c(2, 10), guides = "collect") & theme(legend.position = 'right')
p
# ggsave(filename = "../Submission/Figure 4.png", plot = p, height = 130, width = 140, units = "mm", dpi = 600)

Figure 4: Heat map of the Random Forest regression results showing the ranking of the importance of terrain attributes (based on % increase in Mean Squared Error) in explaining the spatial variabilty of selected colour and geochemical properties within the agricultural and forested sites. Top panel shows an average ranking for each site and across both sites.

4 Discussion

4.1 Variability of soil properties

Variability in soil geochemical properties have been studied at a range of scales including continental (Drew et al., 2010), regional (Rattenbury et al., 2018), watershed (Nanos and Rodríguez Martín, 2012), hillslope/catena, and farm field (Sun et al., 2021b). The objectives of these studies included addressing issues of pollution/contamination, providing benchmark/baseline information, investigating pedological and weathering properties and processes, and soil surveying and mapping (Wilson et al., 2008). Similarly, variability in soil colour, typically using the Munsell colour system, is a commonly reported diagnostic feature used in soil classification and ranges in spatial scale from reconnaissance to detailed soil surveys and maps. For sediment fingerprinting studies, these types of studies are often too site-specific or focus on a smaller subset of soil properties to effectively guide sample design to ensure the desired confidence is met characterizing sources of sediment.

Data distributions in soil science commonly exhibit a positively skewed distribution. This is likely due to several factors including that data of this nature are a semi-bounded distribution, with a lower bound of zero and no upper bound. Hot spots of soil processes, local variations in soil forming factors, and soil/land management practices can also lead to more extreme values (e.g., Vidon et al., 2010). In many cases the cumulative effects of these processes, factors, and practices are multiplicative (i.e., interact) and not linearly additive, resulting in a skewed data distribution. Lastly, the distribution of data will also be a product of the scale of observation, number of samples, and sampling design.

Soil colour properties exhibited a near-normal distribution with a low CV which is consistent with claims that soil hue and value (Munsell colour system) have a low CV (Pennock et al., 2008). These data distribution properties are ideal for statistical and environmental modeling as it typically meets the model assumptions with out requiring transformations. For example, in sediment source fingerprinting, soil properties (i.e., fingerprints) are considered more reliable and robust for use in unmixing models when they show large differences between sources and low variability within each source. Additionally, most mixing models assume fingerprint data are normally distributed. (Luna Miño et al., 2024) demonstrated that soil colour coefficients a*, b*, c*, h*, and x provided good discrimination between the agricultural and forested sites, and the low CV and skewness values reported in Table 1 makes these colour properties ideal fingerprints for sediment source apportionment.

The geochemical properties were more variable and skewed as compared to the soil colour properties. For many trace elements, concentrations are strongly correlated with the proportion of fine-grained material (<2 µm), due to its high specific surface area and enhanced chemical reactivity (Horowitz, 1991). However, in this study the sand-size (>63 µm) material was removed prior to analysis to reduce the effects of grain-size on concentration. This likely resulted in lower variability and less extreme concentrations as compared to other studies that focus on bulk soil samples (<2 mm). In particular, the forested site exhibited a greater amount of variability which is likely due to the more complex topography and geomorphic setting. The floodplain within the forested site is likely accumulating shale-rich material derived from the Manitoba Escarpment which is enriched in trace metals (Nicolas and Bamburak, 2011). This creates a zone of high concentrations relative to upland areas Figure 3. The forested site also had a higher and much more variable soil organic matter content (\(\bar{x}\) = 8.5 %, CV = 51.9 %) as compared to the agricultural site (\(\bar{x}\) = 11.6 %, CV = 16.1 %), which similarly to the grain size distribution, can influence the concentration of many major and trace elements (Horowitz, 1991). These results provide evidence that both land use and landscape complexity both play a role in driving soil property variability.

4.2 Spatial distribution

The difference in the number of soil properties and the magnitude of the spatial auto correlation between the two sites can be used in designing an effect sampling campaign. The agricultural site, which has a simpler topography and a higher degree of spatial autocorrelation, the range values can be used to guide the distance between sampling points and a grid-style sampling regime may be an effective approach. In contrast, the forested site, which has a more complex geomorphic setting and a lower degree of spatial autocorrelation, a stratified sampling design may be the better approach. For example, at the forested site the stratas could include near-stream and hillslope environments. In situations where the soil properties of interest are not known or selected a priori (e.g., sediment fingerprinting) the differences in their spatial autocorrelation are difficult to accommodate in the sampling design. A sampling grid with irregular spacing, including spacing less than 100m, would have provided information on the spatial autocorrelation over shorter distances and reduced the uncertainty in the interpolation of soil properties (Lark and Marchant, 2018).

Mapping the soil properties that have a moderate to high spatial dependence can provide information on underlying soil forming processes and properties. At both sites, to some extent, the patterns appear to reflect the topography of the sites suggesting that geomorphic and hydrologic processes and properties are likely driving the observed patterns. Identifying patterns and understanding the underlying process and properties that drive these patterns are important considerations when designing as soil sampling campaign to successfully meet study objectives, including characterizing soil properties of a field site. In a related context, Koiter (Koiter et al., 2013b) discussed the issues surrounding the use of a statistical only approach to selecting fingerprints and that consideration of how fingerprints have developed improves the robustness of the sediment fingerprinting approach. However, local information on the spatial distribution of geochemical and colour properties at field scales (< 1 km²) is often unavailable, and the processes driving these patterns are also not well documented or studied. When such information does exist, it typically focuses on agronomically important properties (e.g., Mzuku et al., 2005) or is used for soil classification (e.g., Group, 1998). These datasets usually include geochemical properties such as nitrogen (N), phosphorus (P), potassium (K), sulfur (S), calcium (Ca), magnesium (Mg), sodium (Na), iron (Fe), aluminum (Al), nitrate (NO₃⁻), carbonate (CO₃²⁻), bicarbonate (HCO₃⁻), chloride (Cl⁻), and sulfate (SO₄²⁻). They may also include colour characteristics, such as Munsell hue, value, and chroma, as well as other soil properties like texture, organic matter content, and pH. The lack of information on the wide range of soil properties means the researchers are relying on other data, most often elevation, for informing sampling designs.

4.3 Terrain attributes and soil properties

Both the correlation analysis and random forest regression identified elevation as the most influential terrain attribute, followed by the SAGA Wetness Index and relative slope position, in explaining most of the observed variation in soil geochemical and colour properties. This is consistent with the findings of Mashalaba (Mashalaba et al., 2020) who also found that similar terrain attributes were important in predicting a range of other soil properties including texture, bulk density, and hydaulic conductivity. These attributes likely emerged as the most important factor in explaining the observed variability as they are strongly linked to a range of geomorphic and hydrologic process and conditions [(Mello et al., 2022); (Libohova et al., 2024)]. For example, in eroded landscapes in the Prairie region of Canada, Ca concentrations have been found to be higher in upper slope positions from erosion and subsquent exposure of high-carbonate subsoil (Papiernik et al., 2005). In contrast, higher Ca concentrations have been noted in lower slope and depressional areas due to higher solubility of many Ca-minerals (e.g., CaCO₃) and the subsequent downslope transport in solution and reduced leaching losses in these accumulation zones. Landscape position can also have a strong influence on pedogenic process; for example, the translocation of Fe and clay down the soil profile is a diagnostic criteria used in classifying soils (Stonehouse and St. Arnaud, 1971). Soil colour also tends to change in a predictable manner in relation to local relief. Tillage and water erosion results in the net loss of darker organic-rich topsoil from upper slope positions resulting in the exposure of the lighter subsoil (Papiernik et al., 2005). Moisture availability is also greater in the lower slope and depressional areas resulting in increased organic matter production resulting in darker organic-rich topsoil as compared to the upper slope positions. There is also evidence that suggests that soil texture varies with elevation and slope position, with coarser material on upper slopes and finer material accumulating in lower positions (Cox et al., 2003; Kokulan et al., 2018). Given the strong correlation of organic matter and texture with soil geochemistry (Horowitz, 1991) and colour (Viscarra Rossel et al., 2009), these properties may also help explain the observed spatial patterns .

The relative importance of terrain attributes in explaining soil property variability differs both among soil properties and between sites. The land use and the overall geomorphic complexity differences between the two study sites are likely interacting with terrain attributes and influencing the patterns of soil properties and modifying the nature of terrain attribute and soil property relationship. This suggests that these relationships observed in this study may not be broadly generalized. Similarly, information on how terrain attributes influence the spatial distribution of many trace elements and soil colour, beyond the Munsell system, at the field scale is limited in the scientific literature. Additional variables including climate and large-scale landscape features will also influence the observed patterns of soil properties. As a result, using terrain attributes to guide soil sampling or interpret spatial patterns of many soil properties remains challenging.

The impact of sampling design at the field scale on the characterization of soil properties can be substantial (Luna Miño et al., 2024), which in turn can affect the interpretation of data, modeling results, and the conclusion drawn. High-quality LiDAR data or digital elevation models (DEMs) are increasingly openly available in many regions and can be used to create detailed terrain attribute maps. By incorporating terrain attributes into the sampling framework, researchers can ensure that key geomorphic and hydrologic gradients are adequately represented. Ultimately, integrating terrain analysis into sediment source fingerprinting is promising not only as a mechanism to improve the quality of source characterization but to also better link source material to downstream sediment.

5 Conclusions

Understanding the spatial variability and distribution of soil geochemical and colour properties at a field-scale is important for agricultural and environmental research, monitoring, modeling, and management practices. This study conducted both univariate and spatial analyses of a suite of soil geochemical and colour properties at two sites with contrasting land uses. The agricultural site, characterized by gently sloping topography, exhibited lower coefficients of variation, approximately normal data distributions, and moderate to strong spatial autocorrelation across most measured properties. In contrast, the forested site featured more geomorphologically complex terrain, with greater variability in soil properties, data distributions that more frequently deviated from normality, and fewer properties exhibiting spatial autocorrelation. Despite these differences, random forest regression consistently identified elevation, the SAGA Wetness Index, and relative slope position as the three most important terrain attributes explaining the observed variability.

These findings underscore the role of topographic controls on many soil property distributions, regardless of land use. However, the strength and direction of the relationship between terrain attributes and soil property results were inconsistent between both site and soil property. While the study was limited to two sites, the approach demonstrates the value of integrating tools like random forest regression with spatial data to better understand soil-landscape relationships. Future research should expand to broader landscapes and incorporate additional biophysical variables to improve generalizability. Overall, this work highlights how terrain-driven spatial patterns can inform more targeted soil sampling, modeling, and land management strategies.

Acknowledgments

Special thanks and recognition for the field and technical support from A. Avila and the Riding Mountain National Park personnel.

Statements and declarations

Funding

This research was supported by the Natural Sciences and Engineering Research Council of Canada Discovery Grant - From source to sink: Investigating the linkages between sources of sediment and downstream water quality in Canadian watersheds - awarded to AJK (RGPIN-2019-05273).

Competing interests

The authors have no competing interests to declare that are relevant to the content of this article.

Data and code availability

Data and source code for analysis and manuscript available on GitHub: https://github.com/alex-koiter/sampling-design-manuscript

Author contributions

M Luna Miño Methodology; Investigation; Data curation; Formal analysis; Writing - Original Draft; Writing - Review & Editing

A Koiter: Conceptualization; Funding acquisition; Methodology; Investigation; Data curation; Formal analysis; Visualization; Writing - Original Draft; Writing – review and editing; Software; Project administration

T Lychuk: Methodology; Formal analysis; Writing - Review & Editing

A Waddel Methodology; Formal analysis; Writing - Review & Editing

A Moulin Methodology; Formal analysis; Writing - Review & Editing

References

Agisoft, 2021. Agisoft metashape: installer.

Beaudette, D.E., Dahlgren, R.A., O’Geen, A.T., 2013. Terrain-Shape Indices for Modeling Soil Moisture Dynamics. Soil Science Society of America Journal 77, 1696–1710. https://doi.org/10.2136/sssaj2013.02.0048

Beck, H.E., Zimmermann, N.E., McVicar, T.R., Vergopolan, N., Berg, A., Wood, E.F., 2018. Present and future Köppen-Geiger climate classification maps at 1-km resolution. Scientific Data 5, 180214. https://doi.org/10.1038/sdata.2018.214

Boudreault, M., Koiter, A.J., Lobb, D.A., Owens, P.N., Liu, K., Benoy, G., Danielescu, S., Li, S., 2018. Using colour, shape and radionuclide sediment fingerprints to identify sources of sediment in an agricultural watershed in atlantic canada. Canadian Water Resources Journal 43, 347–365. https://doi.org/https://doi.org/10.1080/07011784.2018.1451781

Brown, D.J., Clayton, M.K., McSweeney, K., 2004. Potential terrain controls on soil color, texture contrast and grain-size deposition for the original catena landscape in uganda. Geoderma 122, 51–72. https://doi.org/10.1016/j.geoderma.2003.12.004

Canada, N.R., 2024. High Resolution Digital Elevation Model Mosaic (HRDEM Mosaic) - CanElevation Series.

Collins, A.L., Blackwell, M., Boeckx, P., Chivers, C.-A., Emelko, M., Evrard, O., Foster, I., Gellis, A., Gholami, H., Granger, S., Harris, P., Horowitz, A.J., Laceby, J.P., Martinez-Carreras, N., Minella, J., Mol, L., Nosrati, K., Pulley, S., Silins, U., Silva, Y.J. da, Stone, M., Tiecher, T., Upadhayay, H.R., Zhang, Y., 2020. Sediment source fingerprinting: benchmarking recent outputs, remaining challenges and emerging themes. Journal of Soils and Sediments 20, 4160–4193. https://doi.org/10.1007/s11368-020-02755-4

Collins, A.L., Burak, E., Harris, P., Pulley, S., Cardenas, L., Tang, Q., 2019. Field scale temporal and spatial variability of δ13C, δ15N, TC and TN soil properties: Implications for sediment source tracing. Geoderma 333, 108–122. https://doi.org/10.1016/j.geoderma.2018.07.019

Collins, A.L., Walling, D.E., Leeks, G.J.L., 1997. Source type ascription for fluvial suspended sediment based on a quantitative composite fingerprinting technique. Catena 29, 1–27. https://doi.org/https://doi.org/10.1016/S0341-8162(96)00064-1

Conrad, O., Bechtel, B., Bock, M., Dietrich, H., Fischer, E., Gerlitz, L., Wehberg, J., Wichmann, V., Böhner, J., 2015. System for automated geoscientific analyses (SAGA) v. 2.1.4. Geoscientific Model Development 8, 1991–2007. https://doi.org/10.5194/gmd-8-1991-2015

Cox, M.S., Gerard, P.D., Wardlaw, M.C., Abshire, M.J., 2003. Variability of Selected Soil Properties and Their Relationships with Soybean Yield. Soil Science Society of America Journal 67, 1296–1302. https://doi.org/10.2136/sssaj2003.1296

Drew, L.J., Grunsky, E.C., Sutphin, D.M., Woodruff, L.G., 2010. Multivariate analysis of the geochemistry and mineralogy of soils along two continental-scale transects in north america. Science of The Total Environment 409, 218–227. https://doi.org/10.1016/j.scitotenv.2010.08.004

Du, P., Walling, D.E., 2017. Fingerprinting surficial sediment sources: Exploring some potential problems associated with the spatial variability of source material properties. Journal of Environmental Management 194, 4–15. https://doi.org/https://doi.org/10.1016/j.jenvman.2016.05.066

Ehrlich, W.A., Pratt, L.E., Leclaire, F.P., 1958. Reconnaissance soil survey of west-lake map sheet area.

Environment, Canada, C.C., 2024. Canadian Climate Normals.

Esri, 2024. ArcGIS pro. Esri.

Group, S.C.W., 1998. The Canadian System of Soil Classification, 3rd ed. NRC Research Press, Ottawa, ON, Canada.

Haddadchi, A., Hicks, M., Olley, J.M., Singh, S., Srinivasan, M.s., 2019. Grid-based sediment tracing approach to determine sediment sources. Land Degradation and Development 30, 2088–2106. https://doi.org/10.1002/ldr.3407

Hijmans, R.J., 2024. Terra: Spatial data analysis.

Horowitz, A.J., 1991. A primer on sediment-trace element chemistry, 2nd ed. ed. Lewis Publishers, Chelsea, Michigan, USA.

Kariuki, S.K., Zhang, H., Schroder, J.L., Hanks, T., Payton, M., Morris, T., 2009. Spatial variability and soil sampling in a grazed pasture. Communications in Soil Science and Plant Analysis 40, 16741687.

Koiter, A., 2021. Colour analysis r scripts. https://doi.org/10.5281/zenodo.5123327

Koiter, A.J., Owens, P.N., Petticrew, E.L., Lobb, D.A., 2013a. The behavioural characteristics of sediment properties and their implications for sediment fingerprinting as an approach for identifying sediment sources in river basins. Earth-Science Reviews 125, 24–42. https://doi.org/https://doi.org/10.1016/j.earscirev.2013.05.009

Koiter, A.J., Owens, P.N., Petticrew, E.L., Lobb, D.A., 2013b. The behavioural characteristics of sediment properties and their implications for sediment fingerprinting as an approach for identifying sediment sources in river basins. Earth-Science Reviews 125, 24–42. https://doi.org/https://doi.org/10.1016/j.earscirev.2013.05.009

Kokulan, V., Akinremi, O., Moulin, A.P., Kumaragamage, D., 2018. Importance of terrain attributes in relation to the spatial distribution of soil properties at the micro scale: A case study. Canadian Journal of Soil Science 98, 292–305. https://doi.org/10.1139/cjss-2017-0128

Kuhn, Max, 2008. Building predictive models in r using the caret package. Journal of Statistical Software 28, 1–26. https://doi.org/10.18637/jss.v028.i05

Laceby, J.P., Evrard, O., Smith, H.G., Blake, W.H., Olley, J.M., Minella, J.P.G., Owens, P.N., 2017. The challenges and opportunities of addressing particle size effects in sediment source fingerprinting: A review. Earth-Science Reviews 169, 85–103. https://doi.org/https://doi.org/10.1016/j.earscirev.2017.04.009

Lark, R.M., Marchant, B.P., 2018. How should a spatial-coverage sample design for a geostatistical soil survey be supplemented to support estimation of spatial covariance parameters? Geoderma 319, 89–99. https://doi.org/10.1016/j.geoderma.2017.12.022

Liaw, A., Wiener, M., 2002. Classification and regression by randomForest. R News 2, 18–22.

Libohova, Z., Mancini, M., Winzeler, H.E., Read, Q.D., Sun, N., Beaudette, D., Williams, C., Blackstock, J., Silva, S.H.G., Curi, N., Adhikari, K., Ashworth, A., Minai, J.O., Owens, P.R., 2024. Interpreting the spatial distribution of soil properties with a physically-based distributed hydrological model. Geoderma Regional 39, e00863. https://doi.org/10.1016/j.geodrs.2024.e00863

Lima, F.R.D., Pereira, P., Vasques, I.C.F., Silva Junior, E.C., Mancini, M., Oliveira, J.R., Prianti, M.T.A., Windmöller, C.C., Weindorf, D.C., Curi, N., Ribeiro, B.T., Richardson, J., Marques, J.J., Guilherme, L.R.G., 2023. Predictive modeling of total hg background concentration in soils of the amazon rainforest biome with support of proximal sensors and auxiliary variables. Journal of South American Earth Sciences 129, 104510. https://doi.org/10.1016/j.jsames.2023.104510

Luna Miño, M.A., Koiter, A.J., Lobb, D.A., 2024. Effect of sampling design on characterizing surface soil fingerprinting properties. Journal of Soils and Sediments 24, 2180–2198. https://doi.org/10.1007/s11368-024-03805-x

MacKay, G.H., 1970. A quantitative study of geomorphology of the wilson creek watershed, manitoba (MSC thesis). University of Manitoba, Winnipeg, MB.

Mashalaba, L., Galleguillos, M., Seguel, O., Poblete-Olivares, J., 2020. Predicting spatial variability of selected soil properties using digital soil mapping in a rainfed vineyard of central chile. Geoderma Regional 22, e00289. https://doi.org/10.1016/j.geodrs.2020.e00289

Mello, F.A.O., Demattê, J.A.M., Rizzo, R., Mello, D.C. de, Poppiel, R.R., Silvero, N.E.Q., Safanelli, J.L., Bellinaso, H., Bonfatti, B.R., Gomez, A.M.R., Sousa, G.P.B., 2022. Complex hydrological knowledge to support digital soil mapping. Geoderma 409, 115638. https://doi.org/10.1016/j.geoderma.2021.115638

Mzuku, M., Khosla, R., Reich, R., Inman, D., Smith, F., MacDonald, L., 2005. Spatial variability of measured soil properties across site-specific management zones. Soil Science Society of America Journal 69, 1572–1579. https://doi.org/10.2136/sssaj2005.0062

Naimi, B., Hamm, N. a.s., Groen, T.A., Skidmore, A.K., Toxopeus, A.G., 2014. Where is positional uncertainty a problem for species distribution modelling. Ecography 37, 191–203. https://doi.org/10.1111/j.1600-0587.2013.00205.x

Nanos, N., Rodríguez Martín, J.A., 2012. Multiscale analysis of heavy metal contents in soils: Spatial variability in the duero river basin (spain). Geoderma 189-190, 554–562. https://doi.org/10.1016/j.geoderma.2012.06.006

Nicolas, M.P.B., Bamburak, J.D., 2011. Geochemistry and mineralogy of cretaceous shale, southwestern manitoba (parts of NTS 62F, g, j, k, n, 63C): Phase 2 results. Manitoba Innovation, Energy; Mines, Manitoba Geological Survey.

Papiernik, S.K., Lindstrom, M.J., Schumacher, J.A., Farenhorst, A., Stephens, K.D., Schumacher, T.E., Lobb, D.A., 2005. Variation in soil properties and crop yield across an eroded prairie landscape. Journal of Soil and Water Conservation 60, 8.

Pennock, D., Yates, T., Braidek, J., 2008. Soil sampling designs, in: Carter, M.R., Gregorich, E.G. (Eds.),. CRC Press, Boca Raton, FL, USA, pp. 1–14.

Pulley, S., Collins, A.L., Van der Waal, B., 2018. Variability in the mineral magnetic properties of soils and sediments within a single field in the cape fold mountains, south africa: Implications for sediment source tracing. CATENA 163, 172–183. https://doi.org/10.1016/j.catena.2017.12.019

Pulley, S., Foster, I., Collins, A.L., 2017. The impact of catchment source group classification on the accuracy of sediment fingerprinting outputs. Journal of Environmental Management 194, 16–26.

Rattenbury, M., Martin ,Adam, Baisden ,Troy, Turnbull ,Rose, Rogers, K., 2018. Geochemical baseline soil surveys for understanding element and isotope variation across new zealand. New Zealand Journal of Agricultural Research 61, 347–357. https://doi.org/10.1080/00288233.2018.1426616

Ritchie, J.C., Clebsch, E.E.C., Rudolph, W.K., 1970. Distribution of fallout and natural gamma radionuclides in litter, humus and surface mineral soil layers under natural vegetation in the great smoky mountains, north carolina-tennessee. Health Physics 18, 479.

RStudio, 2024. RStudio: Integrated development environment for r.

Starr, J.L., Meisinger, J.J., Parkin, T.B., 1995. Influence of sample size on chemical and physical soil measurements. Soil Science Society of America Journal 59, 713–719. https://doi.org/10.2136/sssaj1995.03615995005900030012x

Stonehouse, H.B., St. Arnaud, R.J., 1971. Distribution of iron, clay and extractable iron and aluminum in some saskatchewan soils. Canadian Journal of Soil Science 51, 283–292. https://doi.org/10.4141/cjss71-036

Sun, Y., Guo, W., Weindorf, D.C., Sun, F., Deb, S., Cao, G., Neupane, J., Lin, Z., Raihan, A., 2021a. Field-scale spatial variability of soil calcium in a semi-arid region: Implications for soil erosion and site-specific management. PEDOSPHERE 31, 705–714. https://doi.org/10.1016/S1002-0160(21)60019-X

Sun, Y., Guo, W., Weindorf, D.C., Sun, F., Deb, S., Cao, G., Neupane, J., Lin, Z., Raihan, A., 2021b. Field-scale spatial variability of soil calcium in a semi-arid region: Implications for soil erosion and site-specific management. PEDOSPHERE 31, 705–714. https://doi.org/10.1016/S1002-0160(21)60019-X

Team, R.C., 2024. R: A language and environment for statistical computing.

Umali, B.P., Oliver, D.P., Forrester, S., Chittleborough, D.J., Hutson, J.L., Kookana, R.S., Ostendorf, B., 2012. The effect of terrain and management on the spatial variability of soil properties in an apple orchard. CATENA 93, 38–48. https://doi.org/10.1016/j.catena.2012.01.010

Vasu, D., Singh, S.K., Sahu, N., Tiwary, P., Chandran, P., Duraisami, V.P., Ramamurthy, V., Lalitha, M., Kalaiselvi, B., 2017. Assessment of spatial variability of soil properties using geospatial techniques for farm level nutrient management. Soil and Tillage Research 169, 25–34. https://doi.org/10.1016/j.still.2017.01.006

Vidon, P., Allan, C., Burns, D., Duval, T.P., Gurwick, N., Inamdar, S., Lowrance, R., Okay, J., Scott, D., Sebestyen, S., 2010. Hot spots and hot moments in riparian zones: Potential for improved water quality management. Journal of the American Water Resources Association 46, 278–298. https://doi.org/10.1111/j.1752-1688.2010.00420.x

Viscarra Rossel, R.A., Cattle, S.R., Ortega, A., Fouad, Y., 2009. In situ measurements of soil colour, mineral composition and clay content by visNIR spectroscopy. Geoderma 150, 253–266. https://doi.org/10.1016/j.geoderma.2009.01.025

Viscarra Rossel, R.A., Minasny, B., Roudier, P., McBratney, A.B., 2006. Colour space models for soil science. Geoderma 133, 320337.

Wickham, H., 2016. ggplot2: Elegant graphics for data analysis. Springer-Verlag, New York NY U.S.A.

Wilkinson, S.N., Olley, J.M., Furuichi, T., Burton, J., Kinsey-Henderson, A.E., 2015. Sediment source tracing with stratified sampling and weightings based on spatial gradients in soil erosion. Journal of Soils and Sediments 15, 2038–2051. https://doi.org/10.1007/s11368-015-1134-2

Wilson, M.A., Burt, R., Indorante, S.J., Jenkins, A.B., Chiaretti, J.V., Ulmer, M.G., Scheyer, J.M., 2008. Geochemistry in the modern soil survey program. Environmental Monitoring and Assessment 139, 151–171. https://doi.org/10.1007/s10661-007-9822-z

Zhang, S., Huang, Y., Shen, C., Ye, H., Du, Y., 2012. Spatial prediction of soil organic matter using terrain indices and categorical variables as auxiliary information. Geoderma, Entering the digital era: Special issue of pedometrics 2009, beijing 171-172, 35–43. https://doi.org/10.1016/j.geoderma.2011.07.012

Supplemental figures

Figure S1: Summary statistics of all measured colour soil properties at both sites. Error bars represent 1SD and the numeric values indicate the CV.

Figure S2: Summary statistics of all measured geochemical soil properties at both sites. Error bars represent 1SD and the numeric values indicate the CV.

Supplemental tables

Table S1: Description of spectral reflectance colour coefficients used as fingerprints. Reproduced from Boudreault et al. (2018)

Colour space model	Parameter	Abbreviation
RGB	Red	R
RGB	Green	G
RGB	Blue	B
CIE xyY	Chromatic Coordinate x	x
CIE xyY	Chromatic Coordinate y	y
CIE xyY	Brightness	Y
CIE XYZ	Virtual component X	X
CIE XYZ	Virtual component Z	Z
CIE LAB	Metric lightness function	L
CIE LAB	Chromatic coordinate opponent red–green scales	a*
CIE LAB	Chromatic coordinate opponent red–green scales	b*
CIE LUV	Chromatic coordinate opponent blue–yellow scales	u*
CIE LUV	Chromatic oordinate opponent red–green scales	v*
CIE LCH	CIE hue	c*
CIE LCH	CIE chroma	h*

Table S2: Terrain attribute descriptions

Terrain Attribute	Description
Elevation	Meters above sea level
Plan Curvature	Across slope curvature
Profile Curvature	Down slope curvature
SAGA Wetness Index	Similar to the ‘Topographic Wetness Index’ (TWI), but it is based on a modified catchment area calculation, which does not think of the flow as very thin film. As result it predicts for cells situated in valley floors with a small vertical distance to a channel a more realistic, higher potential soil moisture compared to the standard TWI calculation
Catchment Area	Area of upslope contributing area
Relative Slope Position	A value between 0 and 1 illustrating the position of a pixel within the landscape with values approaching 0 indicating streams to pits, and values approaching 1 indicating upper slope positions to peaks
Vertical Distance to Channel	The vertical distance to a channel network base level. The algorithm consists of two major steps: 1. Interpolation of a channel network base level elevation 2. Subtraction of this base level from the original elevations

Table S3: Pearson’s correlation coefficients for soil properties and terrain attributes using interpolated values (10m resolution).

In [13]:

readRDS(file = "./images/correlation.rds")

Property	Elevation	SAGA Wetness Index	Rel. Slope Position	Vert. Dist. Channel	Catchment Area	Plan Curvature	Profile Curvature
Agriculture
Ca	-0.76***	0.59***	-0.26***	-0.25***	0.1***	NS	NS
Co	-0.63***	0.61***	-0.23***	-0.24***	0.14***	NS	NS
Cs	0.58***	-0.45***	0.06***	0.14***	-0.11***	NS	NS
Fe	-0.14***	0.32***	-0.09***	-0.09***	0.09***	NS	NS
Li	-0.4***	0.27***	-0.42***	-0.2***	0.13***	NS	NS
La	0.53***	-0.3***	0.07***	0.15***	NS	NS	NS
Nb	0.48***	-0.52***	0.25***	0.15***	-0.08***	NS	NS
Ni	-0.75***	0.71***	-0.32***	-0.32***	0.2***	NS	NS
Rb	0.81***	-0.71***	0.24***	0.3***	-0.15***	NS	NS
Sr	-0.81***	0.64***	-0.35***	-0.27***	0.12***	NS	NS
a*	0.61***	-0.44***	0.35***	0.23***	-0.12***	NS	NS
b*	0.39***	-0.22***	0.22***	0.1***	-0.09***	NS	NS
c*	0.41***	-0.24***	0.22***	0.11***	-0.09***	NS	NS
h*	-0.22***	0.31***	-0.15***	-0.15***	NS	NS	NS
x	-0.15***	0.25***	-0.23***	-0.12***	NS	NS	NS
Forest
Ca	-0.11***	-0.47***	0.09***	0.4***	0.24***	0.1***	NS
Co	-0.04*	-0.34***	0.05***	0.32***	0.15***	0.05**	-0.03*
Cs	0.09***	-0.4***	0.17***	0.34***	0.18***	0.06***	-0.03*
Li	0.1***	-0.15***	0.05***	0.18***	0.07***	NS	NS
La	NS	-0.23***	NS	0.24***	0.11***	0.05***	NS
Nb	0.12***	0.46***	-0.11***	-0.34***	-0.21***	-0.09***	0.05***
Ni	0.19***	-0.27***	0.14***	0.2***	0.04**	NS	NS
Sr	-0.34***	-0.44***	-0.12***	0.36***	0.26***	0.11***	-0.05**
h*	-0.08***	-0.36***	NS	0.35***	0.22***	0.08***	-0.03*
x	-0.04**	-0.38***	0.09***	0.35***	0.22***	0.09***	-0.03*
*** p < 0.001; ** p < 0.01; * p < 0.05; NS = non-significant at p = 0.05