## 2 Materials and Methods
One of the core methods of the analysis was the machine learning algorithm Random Forest (RF), for which the spatialRF package in R was used [50]. This package is based on the ranger RF package [51]. RF is a popular algorithm since it can find linear and non-linear relationships and patterns in large datasets that contain many features. RF is a large collection of decision trees, individually trained on a random subset of the full training dataset, and tested against the remaining data. This process is called out-of-bag training (OOB) and makes the algorithm robust against overfitting [52]. RF has been used in numerous scientific studies on soil research [53; 54] and metal concentrations specifically [55; 56; 57; 58]. The target variables and features were tested for spatial autocorrelation in the present study. Autocorrelation can lead to errors \begin{table}
\begin{tabular}{c c c c c} \hline \hline
**Category** & **Variable (Unit)** & **Min** & **Max** & **Mean** \\ \hline \multirow{3}{*}{Sacramento River} & River density & 0 & 0.21 & 0.04 \\  & River distance (m) & 0 & 45,137.02 & 9465.59 \\ \hline \multirow{8}{*}{Soil} & Bulk density (g/cm\({}^{3}\)) & 0.98 & 1.72 & 1.32 \\  & AWC (mm/mm) & 0.09 & 0.35 & 0.18 \\  & K (mm/h) & 0.13 & 290 & 13.54 \\  & Carbon cont. (\%) & 0 & 9.88 & 2.09 \\  & Clay cont. (\%) & 10 & 65 & 34.13 \\  & Silt cont. (\%) & 18.02 & 67.71 & 37.71 \\  & Sand cont. (\%) & 9.42 & 67.46 & 28.28 \\  & Rock cont. (\%) & 0 & 22.72 & 0.72 \\  & Soil albedo & 0.01 & 0.23 & 0.03 \\  & USLE K & 0.1 & 0.49 & 0.29 \\ \hline \multirow{4}{*}{Input streams (proximity)} & Stream east 1 (m) & 96,021 & 180,608 & 147,041.9 \\  & Stream east 2 (m) & 98,659.79 & 205,235 & 162,519.3 \\ \cline{1-1}  & Stream east 3 (m) & 89,314.97 & 215,096 & 164,437.1 \\ \cline{1-1}  & Stream west 1 (m) & 77,141.8 & 192,200 & 143,473.7 \\ \cline{1-1}  & Stream west 2 (m) & 134,473.6 & 202,403 & 173,751.1 \\ \hline \multirow{4}{*}{Terrain} & Elevation (m) & \(-\)9 & 44 & 7.12 \\ \cline{1-1}  & Slope (\%) & 0 & 16.57 & 2.47 \\ \cline{1-1}  & TWI & 9.99 & 16.64 & 11.81 \\ \hline \multirow{2}{*}{Vegetation} & EVI & \(-\)0.09 & 0.8 & 0.2 \\ \cline{1-1}  & NDVI & \(-\)0.31 & 0.85 & 0.33 \\ \hline \hline \end{tabular}  
* AWC = available water capacity, K = hydraulic conductivity, USLE K = Universal Soil Loss Equation erosion factor K, TWI = Terrain Wetness index.  
\end{table} 
Table 1: The 22 environmental features used in the respective RF models. in the prediction, and spatialRF incorporates spatial predictors into the model to overcome these issues. Figure 2 gives an overview of the workflow of this study. There were two phases of complexity reduction in the methodology. The first step was to create a model with each of the 38 metal concentrations as the dependent variable and all other metal concentrations as independent features to identify relationships between the distributions of metals. All models that achieved an OOB R2 of >0.7 were selected for further study. The most important variable, i.e., the metal with the closest relationship to the target metal of the respective model, was recorded. This step resulted in 13 individual metal distribution patterns (IMDP). The 13 IMDP were then treated with a Principal Component Analysis (PCA) to reduce complexity further and describe multi-metal distribution patterns (MMDP). The variability of the MMDP was then explained with another set of RF models. In these models, 22 environmental features were included as independent variables. The importance of the variables in these models was used to rank the factors affecting the MMDP.