Remote Sensing of Lake Sediment Core Particle Size Using Hyperspectral Image Analysis

: Hyperspectral imaging has recently emerged in the geosciences as a technology that provides rapid, accurate, and high-resolution information from lake sediment cores. Here we introduce a new methodology to infer particle size distribution, an insightful proxy that tracks past changes in aquatic ecosystems and their catchments, from laboratory hyperspectral images of lake sediment cores. The proposed methodology includes data preparation, spectral preprocessing and transformation, variable selection, and model ﬁtting. We evaluated random forest regression and other commonly used statistical methods to ﬁnd the best model for particle size determination. We tested the performance of combinations of spectral transformation techniques, including absorbance, continuum removal, and ﬁrst and second derivatives of the reﬂectance and absorbance, along with di ﬀ erent regression models including partial least squares, multiple linear regression, principal component regression, and support vector regression, and evaluated the resulting root mean square error (RMSE), R-squared, and mean relative error (MRE). Our results show that a random forest regression model built on spectra absorbance signiﬁcantly outperforms all other models. The new workﬂow demonstrated herein represents a much-improved method for generating inferences from hyperspectral imagery, which opens many new opportunities for advancing the study of sediment archives.


Introduction
Hyperspectral imaging (imaging spectroscopy) has gained attention in the geoscience community because of its detailed representation of sediment features [1]. This technique can be used to determine the composition of sediment cores based on the analysis of absorptions at specific wavelengths of electromagnetic radiation [2]. Due to the rapid advance of scanning technologies and associated algorithms, hyperspectral imaging has emerged in recent years as a tool in paleolimnology, which is focused on tracking changes in a given lake as well as its watershed and airshed [3]. Hyperspectral imaging has shown great promise for the analysis of sediment cores with numerous possible applications [1,4,5], given that it is fast and non-destructive, and represents a high spatial resolution technique (i.e., submillimeter resolution). Butz et al. [4] modified the interpretation of relative absorption band depth (RABD) features, first introduced by Rein and Sirocko [6] for reflectance spectroscopy, and used them to infer algal chloropigments as well as bacterial pigments. In another study, Schneider et al. [7] developed a model for the calibration of sediment core chromatography to

Sample Preparation and Particle Size Measurement
The sediment cores were split lengthwise, with one half used for particle size and other destructive laboratory analyses; the other half was used for hyperspectral imaging and then

Sample Preparation and Particle Size Measurement
The sediment cores were split lengthwise, with one half used for particle size and other destructive laboratory analyses; the other half was used for hyperspectral imaging and then preserved as an archive. Cores used for destructive analyses were sliced into discrete 0.5 or 1.0 cm sections and stored in sealed WhirlPak bags. Sample preparation and analysis of the size distribution of sediments was carried out with laser granulometry using a Horiba granulometer (model LA950v2) following the methods of Narancic et al. [24]. This device operates according to the Mie optical principle (including Fraunhofer approximation) for fine sediments whose size spectrum ranges from 0.11 to 3000 µm. The operation lasted approximately two minutes for each sample, from the introduction of the samples into the granulometer to the export of the results; mean particle size values were calculated based on the Folk and Ward [25] method using the Gradistat program [26].

Hyperspectral Data Acquisition
Prior to image analysis, the acquisition of high-quality sediment images is required, as low quality or substandard images will create significant problems during the subsequent analyses. Sediment core surfaces were cleaned and flattened by removing protrusions (using a knife, single-edge razor blade or a glass microscope slide) in order to get the optimum focus along the detector arrays and expose fine sediment structures, while removing irregularities such as depressions or bumps [27]. Because other planned analyses included light-and temperature-sensitive components such as pigments and DNA, leaving the cores to dry for hours was not possible, and so cores were covered with transparent plastic film to mitigate the effects of water reflectance. Each sediment core was positioned in the middle of the scanning table and in the field of view of the cameras, with the top of the core facing toward the focus grid and the camera, and the level of the core surface was adjusted to the same height as the focus grid. The scanning system uses two rows of seven 50-W quartz-halogen lamps to ensure optimum illumination angle and light intensity, with the light radiated to the surface of the sediment core in all directions (i.e., hemispherical incident radiance) and conical beam reflected back to the sensors [28].
After optimizing the operational parameters such as camera settings and scanning table speed, sediment cores were scanned with a scanner system equipped with two hyperspectral cameras (Specim Ltd., Oulu, Finland) working in the SWIR and VNIR wavelengths of the electromagnetic spectrum. The optical characteristics of the cameras are presented in Table 2. After capturing the hyperspectral data, the normalization of the raw hyperspectral image to the dark and white references was performed according to Equation (1): where HSI Raw is the Raw hyperspectral image, DF av is the dark reference, and WF av is the white reference averaged into a single frame. The output is the fraction of reflected light in relation to the white standard. The analyses were performed on a selected transect in the middle of each sediment core. It is generally assumed that sediment is deposited relatively evenly on lake bottoms at the scale of sediment cores (i.e., typically with diameters <10 cm) and, therefore, that sediment layers are of approximately even thickness. Where present, visually evident laminations in our samples indicated that this was the case. We thus considered that core properties should vary little in cross-section and selected optimum lengthwise transects for two reasons. First, sediment surfaces are not homogenous and usually contain porosities and irregularities; therefore, undisturbed transects result in better model analyses. Second, computational time and memory usage are problematic issues which must be confronted (for 1 m of full sediment core, the data storage capacity exceeds 40 GB of memory).

Noisy Band Removal
Imaging spectroscopy is often corrupted by various types of noise, due mainly to internal sensor malfunction [29]. Hyperspectral sensors are usually affected by noise and have non-linear responses in different spectral bands [30]; one common source is noisy bands with low signal-to-noise ratios (SNR) [31]. Removing the bands that are susceptible to interference by unwanted noise is a crucial step to establishing a robust model between hyperspectral images and other proxies [4]. In this study, we characterized noisy bands by calculating the correlation coefficients between adjacent bands and removing those bands with correlation coefficients below 0.95 [30]. As such, 60 bands of the VNIR image, ranging from 395.82 to 440.58 nm, and eight bands of the SWIR image, ranging from 2538.16 to 2577.33 nm, were recognized as noisy bands and removed.

Spectral Data Denoising
Noisy spectral data can also be generated by external sources, such as high reflectance areas due to water (i.e., saturated regions) or weakly reflected regions because of high porosity, as well as by internal defects inside the instrument such as detector sensitivity or the heat produced by the sensor while it is operational [32]. In order to remove the interference of noise and highlight spectral features, several spectral preprocessing procedures have been proposed, including the standard normal variate [32], the detrend algorithm [33], and the autoscaling method [34]. To remove spectral noise and preserve spectral information, we adopted a three-degree polynomial Savitzky-Golay (SG) smoothing algorithm with window size of 11 points [35], following the recommendations of Vidal and Amigo [36] and Ru et al. [18].

VNIR/SWIR Data Integration
Data fusion techniques have been extensively employed on remote sensing images to fuse data from different sensors to improve model performance [37]. Recent studies showed that combining laboratory SWIR and VNIR images into one integrated image can provide complementary information and increase model accuracy [18,38]. In this study, the fusion of VNIR and SWIR images was performed using the scale-invariant feature transform (SIFT) flow algorithm introduced by Liu et al. [39]. In this Remote Sens. 2020, 12, 3850 6 of 17 procedure, the VNIR image was down-sampled to 200 µm resolution to give it the same pixel size as the SWIR image. Subsequently, the registration was performed by subjecting the first principal components (PCs) of the VNIR and SWIR sensor images to the SIFT flow algorithm [40]. Finally, after transformation of the SWIR image to the VNIR image space and removing the overlapping bands between 969.29 and 1001.49 nm, an integrated hyperspectral image consisting of 932 bands was constructed.

Patch-Wise Calibration Sample Creation
To construct a model relating particle size assays to high-dimensional spectra from sampling image regions containing thousands of pixels (i.e., regions of interest; ROI), the patch-wise approach was used [41]. The purpose of the patch-wise strategy is to build different distributions of the spectral data and to make a high-population training set while incorporating both spectral and spatial information. In this procedure, ROIs were divided into patches of predetermined size and the closest spectrum to the mean spectra was selected as a representative spectrum for each patch. To compute the distance of spectra from the mean spectrum, the Euclidean distance of the cumulative spectrum, which was found to be the most suitable distance function for hyperspectral data, was used [42]. During the calibration process, each representative spectrum was associated with the equivalent measured particle size and treated as an independent sample ( Figure 1, Step 4).

Variable Selection Using Random Forest
Random Forest (RF), first introduced by Breiman [43], is a variation of ensemble classification and regression trees, which consists of many decision trees each built over a random selection of calibration samples from the dataset (bootstrap aggregation) and a random subset of candidate variables [44]. With different tree structures and randomly selected variables, not every tree sees all the features or all the calibration samples, and this guarantees that the trees are not correlated and therefore less prone to over-fitting [45]. During the development of ensemble trees, two-thirds of the bootstrap samples, on average, are used to grow the bagged trees and the remaining samples (out-of-bag; OOB samples) are employed to independently cross-validate the model (see schematic in Figure 2).
A specific characteristic of RF is variable importance (VI), which permits not only evaluation and ranking of variables (i.e., wavelengths) in terms of relative significance but also interpretations of data and identification of the most prominent features [46,47]. There is a vast literature on variable selection using RF and it is still the subject of ongoing research [48,49]. A conservative strategy for RF variable selection is based on the changes of prediction accuracy in a recursive forward addition or backward elimination of variables until the best performance is achieved [48]. These approaches usually get acceptable results but require long computational times [49]. Selection can also be made by advanced but time-consuming approaches such as those of Genuer et al. [50], who introduced a two-step ascending variable selection procedure (i.e., VSURF) for accurate data interpretation and model prediction purposes. Here, to get a stable VI ranking, the RF model was run on all available samples for 25 iterations and the mean scores of importance (i.e., the average importance values of 25 runs) were considered as VI values. The variables of high importance (scores above the average of the VIs) were then selected for data interpretation and model development (see Figure 3).

7
the features or all the calibration samples, and this guarantees that the trees are not correlated and therefore less prone to over-fitting [45]. During the development of ensemble trees, two-thirds of the bootstrap samples, on average, are used to grow the bagged trees and the remaining samples (out-of-bag; OOB samples) are employed to independently cross-validate the model (see schematic in Figure 2). A specific characteristic of RF is variable importance (VI), which permits not only evaluation and ranking of variables (i.e., wavelengths) in terms of relative significance but also interpretations of data and identification of the most prominent features [46,47]. There is a vast literature on variable selection using RF and it is still the subject of ongoing research [48,49]. A conservative strategy for RF variable selection is based on the changes of prediction accuracy in a recursive forward addition or backward elimination of variables until the best performance is achieved [48]. These approaches usually get acceptable results but require long computational times [49]. Selection can also be made by advanced but time-consuming approaches such as those of Genuer et al. [50], who introduced a two-step ascending variable selection procedure (i.e., VSURF) for accurate data interpretation and model prediction purposes. Here, to get a stable VI ranking, the RF model was run on all available samples for 25 iterations and the mean scores of importance (i.e., the average importance values of 25 runs) were considered as VI values. The variables of high importance (scores above the average of the VIs) were then selected for data interpretation and model development (see Figure 3).

Hyperspectral Model Fitting
Prior to fitting the models, we applied different transformation techniques to the raw spectral data to evaluate their performance in model development, including absorbance (AB), continuum removal (CR), reflectance first derivative (RFD), reflectance second derivative (RSD), first derivative of absorbance (AFD), and second derivative of absorbance (ASD) transformations [51,52]. In the training process, the available input spectra were randomly divided into calibration (80%) and

Hyperspectral Model Fitting
Prior to fitting the models, we applied different transformation techniques to the raw spectral data to evaluate their performance in model development, including absorbance (AB), continuum removal (CR), reflectance first derivative (RFD), reflectance second derivative (RSD), first derivative of absorbance (AFD), and second derivative of absorbance (ASD) transformations [51,52]. In the training process, the available input spectra were randomly divided into calibration (80%) and validation (20%) sets, and the calibration set was used for building the predictive models based on RF, SVR, PLSR, MLR, and PCR. The predictive models were generated and validated for the six lakes separately using their respective calibration and validation sets. An aggregate model was also generated based on a combined calibration set which incorporated all six individual datasets to examine the applicability of a single model across different lake environments.
Experiments were conducted using MATLAB and Python programming languages on a 3.4 GHz processor with 16 GB of RAM. To avoid bias, each scenario was repeated 10 times and the averages are reported. All RF regression models were constructed using the parameter values Mtry = n/3, nodesize = 5, and Ntree = 500 (n = number of bands). The grid search method was implemented to find the optimum SVR hyperparameters (kernel function = RBF and C = 1000). The partial least squares regression model was implemented with 10 orthogonal predictors.

Evaluation of the Transformation Techniques
The results of the models developed using each of the datasets showed that the absorbance transformation (AB) was often substantially superior to other transformations in terms of prediction accuracy. In many cases, the other transformation techniques performed more poorly than analyses based on the raw spectra (Table 3). For most of the tested regression models, AB equaled or improved R 2 and RMSE when compared with other transformation techniques. In terms of MRE, AB had the lowest values in the best overall model (using random forest regression; discussed in Section 3.2), but in less performant models other transformation techniques had slightly better MREs.

Evaluation of the Predictive Models
The clear superiority of RF over other regression models is shown by the maximum R 2 values in both calibration and validation sets in all lakes (Figure 4). The RF models were stable and robust in calibration, as R 2 values were generally above 0.90. For the validation set, RF improved R 2 by approximately 21, 19, 11, 34, 9, and 24% over the second-best regression model (SVR) for WIL, JOS, TRU, BEC, STA, and Fury datasets, respectively (Figure 4b).

Evaluation of the Predictive Models
The clear superiority of RF over other regression models is shown by the maximum R 2 values in both calibration and validation sets in all lakes (Figure 4). The RF models were stable and robust in calibration, as R 2 values were generally above 0.90. For the validation set, RF improved R 2 by approximately 21, 19, 11, 34, 9, and 24% over the second-best regression model (SVR) for WIL, JOS, TRU, BEC, STA, and Fury datasets, respectively (Figure 4b).
The high R 2 values that resulted from applying RF models to sediment core reconstructions show the great potential of these models to reproduce distinct particle size variations ( Figure 5). Here, the R 2 values for each dataset are calculated between the average hyperspectral predictions of sampling areas and the measured MPS to validate the inferences of volume information of sediment cores from surface hyperspectral imaging. The comparison of predicted particle size time series with their corresponding measured values clearly shows strong coherence ( Figure 5).  The high R 2 values that resulted from applying RF models to sediment core reconstructions show the great potential of these models to reproduce distinct particle size variations ( Figure 5). Here, the R 2 values for each dataset are calculated between the average hyperspectral predictions of sampling areas and the measured MPS to validate the inferences of volume information of sediment cores from surface hyperspectral imaging. The comparison of predicted particle size time series with their corresponding measured values clearly shows strong coherence ( Figure 5). Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 19 12

Evaluation of the Variable Selection Method
Pairwise comparisons of R 2 values of the three best regression models (RF, SVR, and PLSR) show a good agreement between the models built using the full spectra and selected wavelengths ( Figure 6). Most points are near the 1:1 line in both the calibration and validation scatter plots, indicating similar behavior for the full and selected wavelength models. However, in addition to lower R 2 values, the PLSR and SVR models also had a number of points more divergent from the 1:1 line, suggesting that with these models, stronger performance would be achieved using the full spectral range. By contrast, for RF regression, the performance of the models using selected wavelengths was nearly equivalent to those developed using the full spectral range. The comparison of processing times between the two scenarios (with or without variable selection) clearly shows that using selected wavelengths alone noticeably decreased the computational time, on average by 5, 7 and 8.5 times for SVR, RF, and PLSR, respectively.

Evaluation of the Variable Selection Method
Pairwise comparisons of R 2 values of the three best regression models (RF, SVR, and PLSR) show a good agreement between the models built using the full spectra and selected wavelengths ( Figure 6). Most points are near the 1:1 line in both the calibration and validation scatter plots, indicating similar behavior for the full and selected wavelength models. However, in addition to lower R 2 values, the PLSR and SVR models also had a number of points more divergent from the 1:1 line, suggesting that with these models, stronger performance would be achieved using the full spectral range. By contrast, for RF regression, the performance of the models using selected wavelengths was nearly equivalent to those developed using the full spectral range. The comparison of processing times between the two scenarios (with or without variable selection) clearly shows that using selected wavelengths alone noticeably decreased the computational time, on average by 5, 7 and 8.5 times for SVR, RF, and PLSR, respectively.

Evaluation of the Aggregate Model
It is widely acknowledged that models based on individual datasets more accurately characterize material properties [53]. To be broadly applicable, however, a model should perform well in a variety of different systems. Therefore, after demonstrating the efficacy of the methodology we constructed an aggregate model-integrating all available particle size samples rather than those of individual coresin order to test the performance and applicability for MPS prediction across a variety of lake systems ( Figure 7, Table 4). Individual models outperformed the aggregate model in their respective lakes, but the aggregate model was still reasonably accurate, with an average R 2 of 0.81 across all cores (Table 4). Although the reported R 2 showed an average decline of 0.07 units for the aggregate model, we found that the predictive strength of the aggregate model was stable among the different datasets. While the applicability of the aggregate model may be related to different particle sizes (discussed in Section 4), high R 2 and low RMSE and MRE values confirm that aggregate model is highly reliable (Figure 7).

Evaluation of the Aggregate Model
It is widely acknowledged that models based on individual datasets more accurately characterize material properties [53]. To be broadly applicable, however, a model should perform well in a variety of different systems. Therefore, after demonstrating the efficacy of the methodology we constructed an aggregate model-integrating all available particle size samples rather than those of individual cores-in order to test the performance and applicability for MPS prediction across a variety of lake systems (Figure 7, Table 4). Individual models outperformed the aggregate model in their respective lakes, but the aggregate model was still reasonably accurate, with an average R 2 of 0.81 across all cores (Table 4). Although the reported R 2 showed an average decline of 0.07 units for the aggregate model, we found that the predictive strength of the aggregate model was stable among the different datasets. While the applicability of the aggregate model may be related to different particle sizes (discussed in Section 4), high R 2 and low RMSE and MRE values confirm that aggregate model is highly reliable (Figure 7).

Discussion
The use of hyperspectral images for quantifying MPS requires advanced multivariate analysis techniques, which differentiate the response of the different particle size features through spectral characteristics [46]. In paleolimnology, most existing studies that have developed models for specific materials have given emphasis to the MLR, PLSR, and SVR methods [9,54,55]. However, with the growth of spectroscopy applications and the diversity of sediment characteristics from different catchments, the need for more stable and accurate predictive models is increasing. Based on the accuracy metrics for the total sets, and the R 2 values of the calibration and validation sets, the RF algorithm produced predictive models with higher accuracy and stability compared to other methods, followed by the SVR and PLSR models (Table 3; Figure 4). Moreover, among different transformation techniques, AB transformation-based models revealed a performance superior or equivalent to those models developed using the raw spectral data. In some cases, the behavior of models such as PCR was not consistent among different transformation techniques, and some negative R 2 values for PCR were reported. As R 2 compares the fit of the model with that of a horizontal straight line, negative values happen only when the model does not follow the trend of the data and so fits worse than a horizontal line. Given the superiority of the models, the AB-RF method was applied to the sediment core hyperspectral images to infer mean particle size, which was then compared with the low-resolution time series measured with granulometry ( Figure 5). The high R 2 values for the six datasets (all ≥ 0.72, and 3 of 6 cores ≥ 0.90) together with the highly similar curves, indicate that our proposed hyperspectral image-based methodology can replicate more time-intensive laboratory analyses with high reliability.
Hyperspectral imaging of sediment cores captures information in hundreds of bands, providing the great advantage of building more accurate models using only the most relevant data. Variable importance metrics and the related selected bands (Figure 3) are computed based on the contribution of each band to model development, and they enable the identification of spectral regions that have the greatest impact on the calibration model. In our case, the positions of these key regions are related to the spectral responses of different materials in specific parts of the electromagnetic spectrum [56]. The broad selected variables at wavelengths below 1000 nm may result from iron oxides, carotenoids, and chlorophyll pigments [15,57]. Narrow, well-defined selected bands near 1400 and 1900 nm are mainly associated with the O-H bonds of residual water in organic matter, and the features near 2200 nm may arise from overtones of vibrations related to organic matter or clay minerals [58,59]. However, lake sediments are complex mixtures of different organic and inorganic components, and although the variable selection method efficiently detects those wavelengths that are most strongly associated with changes in MPS, it lacks the capacity to determine the exact characteristics underlying these variations. Comparing different models constructed using all hyperspectral image bands to those with wavelengths selected by the variable importance metric, there was little evident variation in model performance. For SVR and PLSR regression models, accuracies were higher when all bands were used for model fitting because of the greater number of features contributing to the fitting process. However, the results of the RF model using selected wavelengths and full spectra showed almost equivalent accuracies. Therefore, we conclude that the RF variable selection method maintains model accuracy, in addition to the advantages of enabling interpretation and lowering computational time requirements.
The stability of the aggregate model across datasets can be seen for different particle sizes, although the aggregate model appears to underestimate particle sizes above~40 µm, likely due to the limited number of training samples in this size range (Figure 7a,b). As an example, Figure 7c shows the variation of the prediction error with respect to predicted particle size for the Fury dataset, which contains particle sizes above~40 µm. The prediction error increases in regions of the image with larger particle sizes. Overall, the concentrations of the points distributed along the 1:1 line and the values of R 2 , RMSE, and MRE indicate the generally strong concordance between predicted and measured particle size. However, prediction accuracies of the model at particle sizes greater than~40 µm may require further investigation, including assessments of the stability of the RF model following the inclusion of samples that span a wider range of particle sizes. At the other end of the size range, cores with very low values and limited particle size variations, such as BEC (overall range 18-25 µm), may return lower coefficients of determination because prediction errors represent a proportionally larger part of the overall predicted variance. However, even with lower R 2 , the models remain highly accurate when evaluated by RMSE and MRE. We also recommend that the RF model be compared with other particle size measurement techniques, which tend to be designed to make measurements only within a specific range of particle sizes. For example, sieving is a technique where particles greater than 63 µm diameter are commonly lumped together, whereas pipette analysis is only suitable for the analysis of particles less than 63 µm diameter [14]. 15 following the inclusion of samples that span a wider range of particle sizes. At the other end of the size range, cores with very low values and limited particle size variations, such as BEC (overall range 18-25 µm), may return lower coefficients of determination because prediction errors represent a proportionally larger part of the overall predicted variance. However, even with lower R 2 , the models remain highly accurate when evaluated by RMSE and MRE. We also recommend that the RF model be compared with other particle size measurement techniques, which tend to be designed to make measurements only within a specific range of particle sizes. For example, sieving is a technique where particles greater than 63 µm diameter are commonly lumped together, whereas pipette analysis is only suitable for the analysis of particles less than 63 µm diameter [14]. Figure 7. Ensemble scatter plots of (a) measured versus predicted particle size, and (b) residuals of the predicted particle sizes. (c) Graphical representation of predicted mean particle size (MPS) for Fury dataset and equivalent prediction standard error map (middle diagrams; orange line: highresolution MPS profile; blue line: one-centimeter averages).

Conclusions
In this paper, we developed a hyperspectral image-based methodology to deal with the complexity of MPS determination from the information in the VNIR and SWIR portions of the electromagnetic spectrum. We evaluated numerous combinations of transformation techniques and regression models. Based on our results we conclude that: Figure 7. Ensemble scatter plots of (a) measured versus predicted particle size, and (b) residuals of the predicted particle sizes. (c) Graphical representation of predicted mean particle size (MPS) for Fury dataset and equivalent prediction standard error map (middle diagrams; orange line: high-resolution MPS profile; blue line: one-centimeter averages).

Conclusions
In this paper, we developed a hyperspectral image-based methodology to deal with the complexity of MPS determination from the information in the VNIR and SWIR portions of the electromagnetic spectrum. We evaluated numerous combinations of transformation techniques and regression models. Based on our results we conclude that: 1. Among the studied spectral variables, the AB transformation technique produced the best results. It is noteworthy that in some cases the raw spectra (without transformation) performed equally well or improved accuracy compared to AB.
2. The experimental results on the six datasets showed that the RF regression model outperformed all other regression models, with significant improvements relative to the second-best model (SVR). The stability and accuracy of the models constructed with different methods followed the order RF>SVR>PLSR>MLR>PCR.
3. RF variable selection using variable importance metrics reduced the computational burden and increased the interpretability of the data presented, at no cost to the strength of the models. Results from regression models with and without variable selection showed highly similar results.
4. Hyperspectral image-inferred mean particle size distributions in six sediment cores closely reproduced values measured with granulometry, confirming the utility of the technique for rapid, high-resolution environmental reconstructions. The performance of the aggregate model established its reliability and applicability to sediment cores from a variety of aquatic environments.