## Materials and Methods
### Statistical Analysis  
To assess the quality of calibration curves, which directly compare the intensity obtained from the spectrum with the reference elemental concentration of the sample measured by XRF, the Pearson's R coefficient (Equation (2)) and the RMSE (Equation (3)) were used:  
\[R=\frac{\sum_{i=1}^{N}\left(\frac{\frac{y_{i}-\overline{y}}{\sigma_{x}}}{ \overline{y}}\right)\sum_{j=1}^{N}\left(\frac{y_{j}-\overline{y}}{\overline{ \sigma_{y}}}\right)}{N-1} \tag{2}\]  
where \(x_{i}\) and \(y_{i}\) are the wavelength and the LIBS intensity values, respectively, \(\overline{x}\), \(\overline{y}\) and \(\sigma_{x}\), \(\sigma_{y}\) are the mean and standard deviation from the distributions, and \(N\) is the number of points; and  
\[RMSE=\sqrt{\frac{\sum_{i=1}^{N}\left(y^{i}_{ref}-y^{i}_{pred}\right)^{2}}{N}} \tag{3}\]  
where \(y^{i}_{ref}\) and \(y^{i}_{pred}\) are, respectively, the reference and the predicted concentration values for the model, and \(N\) is the number of points. The R value measures the degree of correlation between two variables and ranges from \(-1\) (perfect negative correlation) to \(+1\) (perfect positive correlation). The RSME basically measures the difference between predicted and measured values. Thus, the values of R and RSME provide complementary information for model evaluation. Thus, the analytical performance of each calibration method used for the quantification of a given element can be evaluated by the R value of the best fitting line between predicted and nominal values (R = 1 in the case of perfect agreement), the intercept value (ideally equal to zero), and the RMSE of the prediction (equal to zero in case of perfect agreement). Further, to evaluate the validation curve, the MAPE values was calculated (Equation (4)), which measures the average percentage error of the value predicted by the model compared to the reference value:  
\[MAPE=100\times\frac{1}{N}\sum_{i=1}^{N}\left|\frac{y^{i}_{ref}-y^{i}_{pred}}{ y^{i}_{ref}}\right| \tag{4}\]where \(N\) is the total number of samples, and \(y^{i}_{ref}\) and \(y^{i}_{pred}\) are, respectively, the reference and the predicted concentration values used in the model. The MAPE value provides the quantification error expected for tested samples. The LOD values were calculated using the calibration curve according to the following Equation (5):  
\[LOD=\frac{3\times\sigma}{\alpha} \tag{5}\]  
where \(\sigma\) is the error of the linear coefficient, and \(\alpha\) represents the slope of the calibration curve. The LOD represents the minimum concentration value of the analyte that can be measured by the technique used, i.e., the technique sensitivity as the detection of concentrations below the LOD are not possible. In order to avoid possible overfitting, which may occur if data are adjusted very well to the training data but exhibit a low predictive capacity, they were validated by the LOO-CV method by which one sample is taken out from the training set and used later for testing. In this way, a concentration can be predicted for the external sample and compared with the reference value. The process was repeated by leaving each sample out of the training so that all samples could be validated by the model trained with the remaining _N-1_ samples. Although this validation process is more complicated computationally, it can achieve a result more reliable than that obtained by the k-fold validation method, as the former one is independent from the random selection of the training set. As the test sample does not participate in the training set, it works as an external measure. The only relationship between the test sample and the training data set was that the measurements were carried out in a similar way (same day and equipment). Although the ideal procedure would be to validate the protocol using measurements in situ and perform it on different days, this is strictly required only in the final stages of implementation of the equipment for industrial applications.