## 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 [@starr1995; @kariuki2009].

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 [@pulley2017]. However, investigating spatial variability at smaller scales is less common [e.g., @du2017; @pulley2018; @collins2019; @lunamiño2024] and remains a research priority [@collins2020].

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 [@wilkinson2015]. 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 [@lunamiño2024].

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 [@wilkinson2015; @du2017] and using spatially interpolated maps of fingerprint values to provide a finer resolution of the fingerprint variability within each source [@haddadchi2019].

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 [@collins1997]. However, there are concerns that this approach may result in the inclusion of false positives (i.e., type I error) or non-conservative fingerprints [@koiter2013]. 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 [@collins2020]. 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 [@lunamiño2024]. 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., ^137^Cs, ^7^Be; @ritchie1970], plant nutrients [e.g., N, P; @vasu2017], soil colour [e.g., hue, value; @viscarrarossel2006], major non-acid forming cations [e.g., Ca, Na; @sun2021]. 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 [@beaudette2013], texture [@kokulan2018], colour [@brown2004], organic matter [@zhang2012], conductivity [@umali2012], and geochemistry [@lima2023]. 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 [-@lunamiño2024] 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.

## Methods

### 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 (@fig-location_map). 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 [@ehrlich1958]. 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 [@ehrlich1958].

The Köppen-Geiger climate classification of the WCW is cold, without dry season, and with warm summer (Dfb) [@beck2018]. The average annual precipitation is \~539 mm, with approximately 27% falling as snow with a mean annual temperature is 3.0°C [@environmentandclimatechangecanada2024]. The hydrology of the watershed is snowmelt dominated with \~ 80% of the cumulative runoff occurring during the spring season (May and June) [@mackay1970].



### Soil sampling and analysis

This study uses samples and data collected as part of the grid sampling design outlined in Luna Miño [-@lunamiño2024]. 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 [@laceby2017]. 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. [-@boudreault2018], fifteen colour coefficients (R, G, B, x, y, Y, X, Z, L, a\*, b\*, u\*, v\*, c\*, h\*) were calculated for each sample [@koiter2021]. Based on the results of Luna Miño [-@lunamiño2024], 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.

### Geospatial and terrain analysis

All geostatistics were performed with ArcGIS Pro [v 3.3.0 @esri2024]. 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 [@naturalresourcescanada2024]. 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 @agisoft2021]. Ordinary kriging was used to calculate a 1 m gridded digital elevation model for each site. Geographic information software [SAGA v2.1.4 @conrad2015] 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 (@supptab-terrain).

### Data analysis

All subsequent statistical analysis was undertaken using R statistical Software v4.4.0 [@rcoreteam2024] through RStudio Integrated Development Environment v2024.04.2 [@rstudio2024]. Plots and maps were created using the R package `ggplot2` v `r packageVersion("ggplot2")` [@wickham2016]. 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` v`r packageVersion("terra")` @hijmans2024]. Random Forest Regression [`randomForest` v`r packageVersion("randomForest")` @liaw2002] 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` v`r packageVersion("usdm")` @Naimi2014]. 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` v`r packageVersion("caret")` @kuhn2008]. 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^~2~^ 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.

## Results

### 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 (@tbl-univariate-summary). 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 (@tbl-univariate-summary). 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.



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 (@tbl-univariate2-summary). 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.



### 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 (@tbl-geocol-semivariogram). 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 (@fig-ag_map). 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 (@tbl-geocol-semivariogram). 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 (@fig-forest_map).







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 (@supptab-correlation2). 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 (@tbl-rf-summary). 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 (@fig-rf-results). 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.





## Discussion

### Variability of soil properties

Variability in soil geochemical properties have been studied at a range of scales including continental [@drew2010], regional [@rattenbury2018], watershed [@nanos2012], hillslope/catena, and farm field [@sun2021a]. 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 [@wilson2008]. 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. 

<!-- The increase use of soil spectroscopy [@poppiel2022] and creation of soil spectral libraries [ e.g., @canadian] will allow researchers to also investigate soil colour variability using other colour space models [@viscarrarossel2006].  -->

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., @vidon2010]. 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 [@pennock2008]. 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. @lunamiño2024 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 @tbl-univariate-summary 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 [@horowitz1991]. 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 [@nicolas2011]. This creates a zone of high concentrations relative to upland areas [@fig-forest_map]. 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 [@horowitz1991]. These results provide evidence that both land use and landscape complexity both play a role in driving soil property variability.

### 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 [@lark2018].

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 [-@koiter2013a] 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., @mzuku2005] or is used for soil classification [e.g., @soilclassificationworkinggroup1998]. 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.

### 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 [-@mashalaba2020] 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 \[@mello2022; @libohova2024\]. 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 [@papiernik2005]. 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~3~) 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 [@stonehouse1971]. 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 [@papiernik2005]. 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 [@kokulan2018; @cox2003]. Given the strong correlation of organic matter and texture with soil geochemistry [@horowitz1991] and colour [@viscarrarossel2009], 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 [@lunamiño2024], 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.

## 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 {.unnumbered}

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

## Statements and declarations {.unnumbered}

### Funding {.unnumbered}

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 {.unnumbered}

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

### Data and code availability {.unnumbered}

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

### Author contributions {.unnumbered}

**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 {.unnumbered}

::: {#refs}
:::

## Supplemental figures {.unnumbered}

::: {#suppfig-colour_summary}
![](images/colour_summary.png){width="100%"}

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

::: {#suppfig-geo_summary}
![](images/geo_summary.png){width="100%"}

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

## Supplemental tables {.unnumbered}

::: {#supptab-abbrev}
+------------------------+---------------------------------------------------+---------------+
| **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\**         |
+------------------------+---------------------------------------------------+---------------+

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

::: {#supptab-terrain}
+------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
| 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                                                                                                                                                                                                                                                                                                  |
+------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+

: Terrain attribute descriptions
:::

::: {.content-visible unless-format="pdf"}
::: {#supptab-correlation2}
```{r}
readRDS(file = "./images/correlation.rds")
```

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

<!-- ## Random Forest Importance Data -->

<!--  -->

<!--  -->

<!--  -->

<!--  -->

<!--  -->

<!--  -->

<!--  -->

<!--  -->

<!--  -->

<!--  -->

<!--  -->

<!--  -->

<!--  -->

<!--  -->

<!--  -->

<!--  -->

<!--  -->

<!--  -->

<!--  -->

<!--  -->

<!--  -->

<!--  -->
