Establishing a new workflow in the study of core reduction intensity and distribution

: New methodological approaches focused on studying the reduction and use-life of stone tools have emerged in recent years, enabling researchers to move beyond strict technical and technological characterizations and explore specific aspects of occupation dynamics and economic management of resources. Previous studies have shown the importance of reduction distributions of individual measurements rather than averaged values. In this sense, survival analysis, and more specifically Weibull distributions, are one of the main inferential tools used in reduction studies. However, the resolution of Weibull distribution obtained from different methods has not been tested experimentally. In this paper, we present an evaluation of some of the main methods used in the study of core reduction intensity, such as the Volumetric Reconstruction Method, the Scar Density Index, and the non-cortical surface percentage. Our results show 1) strong and positive correlations between these approaches and actual reduction intensity, 2) similar Weibull distributions for non-cortical surface percentage, Volumetric Reconstruction Method, and logarithmic transformation of Scar Density Index. In addition, 3) the results from each method show a similar intra-assemblage variation, with a high percentage of agreement between them. As a result, all the evaluated proposals are useful and reliable methods for estimating the degree of reduction. Finally, a workflow is proposed for approaching reduction in archaeological assemblages by integrating different methods in the same study.

New methodological approaches focused on studying the reduction and use-life of stone tools have emerged in recent years, enabling researchers to move beyond strict technical and technological characterizations and explore specific aspects of occupation dynamics and economic management of resources. Previous studies have shown the importance of reduction distributions of individual measurements rather than averaged values. In this sense, survival analysis, and more specifically Weibull distributions, are one of the main inferential tools used in reduction studies. However, the resolution of Weibull distribution obtained from different methods has not been tested experimentally. In this paper, we present an evaluation of some of the main methods used in the study of core reduction intensity, such as the Volumetric Reconstruction Method, the Scar Density Index, and the non-cortical surface percentage. Our results show 1) strong and positive correlations between these approaches and actual reduction intensity, 2) similar Weibull distributions for non-cortical surface percentage, Volumetric Reconstruction Method, and logarithmic transformation of Scar Density Index. In addition, 3) the results from each method show a similar intra-assemblage variation, with a high percentage of agreement between them.
As a result, all the evaluated proposals are useful and reliable methods for estimating the degree of reduction. Finally, a workflow is proposed for approaching reduction in archaeological assemblages by integrating different methods in the same study.

Introduction and background
Lithic technology, understood from a broader perspective in which the management of tools does not only reflect the technological systems of human groups, but also their territorial and economic organisation, has made it possible to study the relationship between both the management of raw materials and the lithic tools with environmental variables. Thus, through different technological approaches we can delve into behavioural aspects of past human groups such as adaptability to environmental conditions (i.e., the proximity and availability of raw materials), occupational patterns (i.e., the type, intensity and duration of occupations), transport and mobility patterns, among others (e.g., Andrefsky 1994;Bamforth 1986;1990;Blades 2003;Kelly 1992;Kuhn 1991;Marreiros & Bicho 2013;McCall 2012;Morales 2016;Nelson 1991;Parry & Kelly 1987).
Some authors have suggested that the intensity of reduction of any archaeological assemblage strongly depends on the characteristics of the occupation (settlement and mobility patterns, site function, etc.) and the environment (availability and proximity of raw materials, etc.), or is even involved in site formation processes (Schiffer 1987). Thus, depending on the characteristics of the occupation and the environment, there would be different resource management patterns that can be manifested in the raw material provisioning, transport strategies, tool maintenance or re-sharpening, and discard or replacement of stone tools, among others (Bamforth 1986;Binford 1979;Kuhn 1991;Odell 1996;Shott 1986;. The growing interest in these issues explains the great deal of work carried out on the study of reduction intensity in the last years. These studies have looked for variables that correlate with the degree of reduction on different artifact categories, such as cores (Clarkson 2013;Cueva-Temprana et al. 2022;Douglass et al. 2018;Lombao 2021;Lombao et al. 2020;Reeves et al. 2021), handaxes and other bifacial tools (Clarkson 2002;Li et al. 2015;Shipton et al. 2013;Shipton & Clarkson 2015a;, or discrete typologies such as distally retouched end-scrapers (Eren et al. 2005;Morales 2016;Morales et al. 2015a;Shott & Seeman 2015), backed blades (Muller et al. 2018), Aterian tanged points (Ioviţǎ 2011), or unifacial lateral scrapers (Bustos-Pérez & Baena 2019Eren et al. 2005;Eren & Prendergast 2008;Eren & Sampson 2009;Hiscock & Clarkson 2009;Kuhn 1990).
This proliferation of approaches reflects the urgent need to adapt the methodologies to each object of study (e.g., type of tool) or to the addressed problem (e.g., raw material management, mobility patterns, etc.) (Andrefsky 2009). However, depending on their theoretical basis and empirical approach, these methods return results in disparate scales, which in some cases are even expressed as intangible values difficult to correlate intuitively with the percentage of volume extracted. Therefore, assuming the impossibility of attaining a universal reduction index to handle all the technological and typological variability from the archaeological record, efforts should focus on the search for analytical tools to test their effectiveness, to compare the results from different approaches and to build more robust inferences by comparing and integrating the results obtained from different approaches.
Furthermore, some authors (Shott 2002) have emphasised the importance of reduction distributions over averaged values, as they can provide more information in those assemblages whose reduction intensity distributions are asymmetric or non-parametric. One of the main statistical approaches in this regard are the Weibull distributions, popularised by Shott (Shott 1996;2002;2010;Shott & Sillitoe 2005) in lithic studies, and subsequently applied by other authors (Cardillo et al. 2021;Douglass et al. 2018;Lin et al. 2016;Lombao et al. 2022;Morales 2016). The Weibull distribution is a statistical test that models failure in relation to time. In reduction studies, failure can be equated with the moment of discard and time with the reduction intensity, use-life or curation (Shott 1996).
Although the potential of such approaches in reduction studies has been widely explored in archaeological case studies, it should be noted that when it comes to experimentally assessing the validity of a reduction index, regression and determination coefficients have been more widely used (Hiscock & Tabrett 2010;Lombao et al. 2019). However, not much attention has been devoted to contrasting the types of Weibull distributions obtained through different reduction methods and the real reduction distributions of given lithic assemblages, hindering the interpretation of the assemblage's reduction pattern inferred from different methods.
For these reasons, the aim of this paper is to review different methodologies focused on the study of core reduction intensity by comparing Weibull distributions. These methodologies have been measured on a controlled experimental assemblage of reduced cores. This will allow to evaluate the inferential potential of the different methodologies. In addition, the combination of different approaches to a single study case are explored.

Materials
In this study a set of experimental cores produced for previous works (Lombao et al. 2020) has been used The experimental assemblage is composed of 64 riverine cobbles collected from the terraces of the Arlanzón river, comprising a wide variability of morphologies and sizes (381-4424g). The knapping has been carried out by four different knappers, using free-hand percussion with hard hammers in all cases, and following three knapping strategies (unifacial unipolar, bifacial multipolar centripetal and multifacial multipolar) as well as handaxe production ( Figure 1). Each knapping strategy was used on 16 blanks. The procedure was not restrictive by any means, trying to generate variability in the reduction intensity, (i.e., in the percentage of volume or mass lost through reduction). In this way, the 64 blanks have been reduced covering a range of between 10% and 95% reduction. The 3D models of the cobbles prior to knapping and the resulting cores are available at Zenodo (Lombao 2019).

Reduction intensity methods
In this study, we employed the Scar Density Index, Percentage of non-cortical surface, and Volumetric Reconstruction Method to estimate the intensity of reduction. Each method has its unique characteristics, which determine its suitability for archaeological assemblages and the range of values it yields.
1. The Scar Density Index (SDI) is based on the relative increase in the density of flake scars regarding the core's surface area (size) since the latter decreases as reduction advances. Therefore, this index is based on the relation of the number of flake scars larger than 10mm on a core to its surface area (cm 2 ) as a proxy for reduction intensity (Clarkson 2013).
2. The percentage of non-cortical surface (PerNCO) is obtained by dividing the noncortical surface area by the total surface area of each core based on the premise that as reduction progresses the amount of cortical surface decreases (Li et al. 2015). The quantification of cortical and non-cortical surface areas has been measured on the 3D models of the cores, since it allows greater precision (Lin et al. 2010 Lombao et al. 2022), although it has also been explored in other technological contexts such as Mode 2 (Acheulean) or in the Mode 2-3 transition (Acheulean-Mousterian transition) (Lombao 2021). The technical analysis of the cores was used to estimate the number of flaking generations on the core's maximum morphometric axes, to define the number of corrections to be applied to each maximum dimension. When n flaking generations were identified along the maximum length and width, n times the mean platform thickness from the assemblage's flakes was added to the corresponding dimension. For the maximum thickness, the mean flake thickness was used. The median was used instead of the mean when the distribution of the measured flake thickness or platform thickness was non-parametric.

Statistical Procedures
The statistical analysis of the data has been carried out following three steps: 1) exploratory data analysis, 2) correlation analysis, and 3) Weibull distributions.

Exploratory data analysis
First, the results from each method were tested for similarities through the central tendency estimators, such as the mean, median, standard deviation (SD), coefficient of variation (CV), and interquartile range (IQR). Then, this data was contrasted through the Kruskal -Wallis (K-W) method for the variance in the medians in non-parametric distributions.
Individual histograms and density lines were produced to explore the shape of each distribution. Then, the skewness (the degree of symmetry of a distribution) was calculated using the Pearson's coefficient of Skewness (PSI). This index is obtained by multiplying the difference between the mean and median by three, and then dividing the result by the standard deviation Furthermore, the kurtosis was also calculated to measure the shape and the spread of the probability distribution. Finally, the Shapiro-Wilk test (S-W) was used to evaluate the parametric character of distributions and the Kolmogorov -Smirnov (K-S) test applied to test for differences in the distributions.

Correlation analysis
Secondly, a linear correlation and Pearson's R test were carried out between the real reduction data and the reduction estimated by each of the methods. In addition, the correlation between the different methods was also tested.

Weibull distributions
Finally, the two parameters defining the Weibull distribution were used to compare the survivorship curves provided by each proxy. The scale parameter (α) represents the point at which 63.2% of the sample has failed (Dorner 1999), while shape parameter (β) defines the shape of the hazard function of the distribution, so that different values in β produce different slopes on a probability plot (Shott 2002;Shott & Sillitoe 2004).

Data transformation
Considering the exponential correlation documented by Clarkson (2013) between the SDI and the original cores' remaining mass and its logarithmic conversion to obtain a linear correlation (Clarkson 2013: 4351), the SDI data was log-transformed (LogSDI) to explore its effect on the data distribution. As the logarithmic transformation implies obtaining negative values, the values have been converted into positive by multiplying by -1.
To level off the α parameter in the Weibull distributions obtained by the different methods, the values from SDI and LogSDI have been rescaled between 0 and 100 (min-max normalization). Since it is not possible to apply the Weibull distributions to datasets with a value of 0, cores with a standardized value of 0 have been discarded for calculating the Weibull distributions.

Sub-samplings
To test the performance of these transformations in estimating Weibull distributions, three sub-samples from the original experimental dataset were performed to simulate two theoretical reduction scenarios and one archaeological situation.
In the first scenario, we sought to generate a "long-lived" discard pattern (i.e., a discard rate characterised by a slow drop-out in the initial phases), which accelerates as the reduction progresses. For this purpose, all cores with a real reduction percentage (% of extracted volume) higher than 50% have been selected.
The second scenario aims to simulate a situation with a high discard ratio during the initial phases that decreases as reduction advances. For this case, all the cores from the experimental dataset showing a reduction percentage lower than 30% and higher than 90% were retained.
In addition, a third sub-sampling was performed with a two-folded purpose: to evaluate the sensitivity of the Weibull distributions obtained through the different methods to random variations in the core assemblages, and to evaluate how the Weibull distributions change as resolution is lost within the same assemblage. Thus, random sub-samples of the experimental assemblage dataset were performed sequentially forcing a progressive loss of the 10% of the assemblage after each step.

Intra-assemblage variation
Finally, the effectiveness of each method in analysing intra-assemblage reduction variation was measured. Three data clusters were calculated for the real reduction values distribution using the k-means function from the "stats" R package (v. 4.3.0), using the algorithm of Hartigan and Wong (1979). Then, three data clusters were also established for the values produced by each method (VRM, PerNCO and SDI) following the same procedure to calculate the percentage of coincidence between clusters. Finally, a Principal Component Analysis combining the different reduction analysis methods was performed to evaluate the resolution of the multi-method approach in distinguishing the original three groups. PCA analysis and visualization was performed using the "FactoMineR" (v.2.7) (Lê et al. 2008) and "Factoextra"

Reproducibility
All figures and statistical analyses have been performed in R environment (R Core Team 2021; R Studio Team 2020). Following the proposals for transparency in research (Marwick 2017), we include the R code as well as the dataset and supplementary materials in Zenodo ). Table 1 shows that the VRM and PerNCO present similar mean and median to the real percentage of extracted volume, and no statistically significant differences were found (K-W (p) > 0.05). As the SDI and LogSDI provide values on a different scale, a comparison without data transformation is not possible. Regarding the dispersion of the values, the Coefficient of Variation (CV) indicates that both the VRM and PerNCO present similar distributions than the real reduction data, contrasting with the values obtained through the SDI, with a remarkably higher CV, and for the LogSDI, with a lower CV ( Table 1).

Descriptive statistics
The distribution of the percentage of real extracted volume shows a very marked symmetry ( Figure 2), reflected by the low Skewness value. These values are similar to those from the VRM, although in this case there is a certain bimodal character. The PerNCO and the LogSDI show a more pronounced negative skewness, contrasting with the considerable positive skewness shown by the SDI values. Despite these differences all the datasets show parametric distributions (Shapiro-Wilk (p) > 0.05), except for the LogSDI (Table 1). The comparison of the distributions (Kolmogorov-Smirnov) shows statistically significant differences between both the SDI and LogSDI and the real reduction data, something that does not occur for the VRM and the PerNCO (see Supplementary File 1).

Correlation analysis
The correlations between each method and the real reduction data are positive, high, and statistically significant, with the VRM values standing out (Pearson's R= 0.85, r 2 = 0.72, p < 0.05). The PerNCO, SDI and LogSDI deliver slightly lower values (Figure 3), although they are still strong and statistically significant positive correlations (Supplementary File 2). Furthermore, the high correlation between the different methods stands out, with Pearson's Rvalues not lower than 0.73 in any case. This indicates a similar performance despite the differences observed among the distributions.

Weibull Distributions
The results show that the VRM is the method that presents the most similar survival curve to the real data, with remarkably similar values for the Shape and Scale parameters, followed by the PerNCO. In contrast, both SDI and LogSDI show values diverging from the real percentage of extracted volume ( Table 2). The Scale parameter differs in the case of SDI and LogSDI because their values are calculated in reference to different scales than the other methods (ranging from 0 to 100). To overcome this situation, the SDI and LogSDI values were rescaled to 0-100.
This normalization does not only affect the Scale values, but the Shape ones, generating slightly lower values. After the rescaling, LogSDI presents a survivorship curve similar to the real one. The SDI, however, presents significantly lower values for both Shape and Scale, which translates into a discard pattern characterised by a greater discard ratio in the initial phases that decreases as the reduction advances (Figure 4a, b).  4. a) Distribuciones de Weibull para el conjunto experimental completo tras la normalización SDI y LogSDI min-max; b) Gráfico bivariante que muestra la relación entre los valores de shape y scale de las distribuciones de Weibull para el conjunto experimental completo tras la normalización SDI y LogSDI min-max; c) Distribuciones de Weibull para la submuestra 1; d) Distribuciones de Weibull para la submuestra 2; y e) Gráfico bivariante que muestra la relación entre el valor de forma y escala de las distribuciones de Weibull para las dos submuestras ). (2023)  To check whether the data transformation is effective in different scenarios, the parameters of the Weibull distribution were calculated for the two subsamples simulating two drastically different discard patterns. VRM, PerNCO and LogSDI allowed to detect these differences among the two subsamples. However, through the SDI, it is not possible to identify major differences in the survivorship curves obtained for each simulation, resulting in low reduction scenarios (Figure 4c-e).

Journal of Lithic Studies
Finally, to evaluate the sensitivity of the Weibull distribution curves, the experimental assemblage was randomly and sequentially subsampled into assemblages 10% smaller, allowing to observe the evolution of the Shape and Scale parameters from both the real data and the different methods as the assemblage composition is randomly transformed.
The results do not show large variations in the Shape and Scale values of the real reduction data as the percentage of cores in the assemblage decreases. Similarly, the PerNCO and VRM remain stable throughout the simulations. However, the LogSDI min-max transformation shows larger variations of the Shape value, while the Scale value is more stable. In contrast, the SDI is less sensitive to these changes, although it systematically shows lower values for both Shape and Scale. Despite differences on the correlation slopes for each method, there is no statistically significant increasing or decreasing tendency for the Shape and Scale values as the simulated resolution of the assemblage is reduced ( Figure 5, Supplementary File 3).  ). Figura 5. Gráfico bivariante que muestra la evolución de los valores de shape (a) y scale (b) a medida que disminuye el porcentaje de núcleos, y regresiones obtenidas entre los valores de shape (c) y scale (d) a medida que disminuye el porcentaje de núcleos analizados ). (2023)

Intra-assemblage variation
To further evaluate the relationships between each reduction index and the real reduction, three data clusters were calculated from the real data and each reduction index. The results indicate a 70.3% coincidence in cluster assignment for the VRM compared to the clusters from the real data, while coincidence drops to 60% in the case of the PerNCO and to 57.8% in the case of SDI and LogSDI. In all three methods, the highest match rate occurs in the lowest reduction group, while the match percentage decreases progressively in SDI and VRM, reaching the lowest match values in the higher reduction stages (Table 3). In PerNCO, on the other hand, the percentage of coincidence is lower in intermediate stages ( Figure 6).
Despite these slight discrepancies among the clustering, the PCA allows to distinguish the three clusters established based on the real reduction data from the combination of the three methods (Figure 7). In the PCA, Dimension 1 explains 85.6% of the variance, while Dimension 2 explains 8.7%, and Dimension 3 explains 5.7%. This indicates very similar results in each core regardless of method (Dimension 1), as all three proxies used move in the same direction along the x-axis, while Dimension 2 would reflect the discrepancies between LogSDI and the other two methods, as they move in the opposite direction along the y-axis. However, there is some overlap between the extreme groups (lower and higher reduction) and the intermediate group.  ). Figura 6. Gráfico que muestra la agrupación obtenida mediante los diferentes métodos en relación con los datos reales de reducción. A) Datos reales; B) VRM; C) LogSDI; D) Porcentaje de superficie no cortical ). (2023)   Groups 1, 2 and 3 correspond to the clustering obtained from the real reduction data. The ellipses show 95% confidence interval around group mean points ). Figura 7. Gráfico del análisis de componentes principales que muestra la dispersión de los núcleos en los tres métodos analizados. Los grupos 1, 2 y 3 corresponden a la agrupación obtenida a partir de los datos reales de reducción. Las elipses muestran el intervalo de confianza del 95% alrededor de los puntos medios de los grupos ).

Discussion
In this paper, we explored the performance of different methods aimed at quantifying reduction intensity in cores. For the first time, this exploration has been carried from a two-fold perspective: by comparing the resolution of the Weibull distribution curves obtained through each method and by studying the intra-assemblage reduction variability at the individual level using multivariate statistics to combine the different reduction indexes.
The correlations obtained between the different methods, as well as with the real data, indicate that all the approaches evaluated here constitute useful and reliable methods for estimating the degree of core reduction. This is observable as they operate in an apparently similar way, despite the differences in their scales. However, the value distribution affects the resolution of the Weibull distribution, as we found significant differences in the survival curves. Thus, while the PerNCO, and especially the VRM, present values of Shape and Scale very similar to the real reduction data, the SDI presents a lower value of Shape.
Therefore, an inversely proportional relationship can be observed between applicability and performance or, in other words, between generalisation and effectiveness. Thus, VRM and PerNCO performed better in this experiment, providing more similar distributions and higher correlations with the real reduction data. Both methods also present the same time intrinsic requirements that make their generalised application on archaeological assemblages more difficult. For example, VRM performance has not been tested on cores on flake. Moreover, it requires the preliminary evaluation of flakes from the assemblage, limiting its applicability in those sites where there are no associated flakes (Lombao et al. 2020) or where the initial raw material formats are not well known (Cueva-Temprana et al. 2022;Lombao et al. 2022). The PerNCO, on the other hand, also presents limitations on cores on flake, as well as being affected by the knapping strategy, since the loss of cortical surface occurs at a different rates depending on the volumetric evolution of the knapping .
As for the SDI, its simplicity has led to a more widespread use, as demonstrated by its application in many chronologically, geographically, and technologically diverse archaeological contexts ( However, in this study, we have documented quite different results regarding the real distributions. Thus, the disparity in the results obtained on Weibull distributions based on the SDI limits their application to archaeological contexts and suggests that it is necessary to review the way in which the effectiveness of reduction methods is analysed. The evaluation of the effectiveness of the reduction methods has progressed in the recent years (Hiscock & Tabrett 2010;Lombao et al. 2019). This has been reached by estimating the effect of some variables on different indexes and applying metrics common in other fields of study. That is the case of using the Average Error, the Mean Average Error, or the Root Mean Square Error (Lombao et al. 2020); or testing the differences between predicted and real reduction between reduction intensity intervals and original blank size (Douglass et al. 2018;Reeves et al. 2021). These analyses have allowed further evaluation of the behaviour of the indexes along the reduction process to address the effectiveness of continuous variables.
The application of new analytical methods has allowed a more detailed assessment of the qualities of an ideal reduction index (Hiscock & Tabrett 2010). The analysis of correlations and especially the coefficient of determination (r-squared) allow to evaluate the inferential power of the reduction indexes, while the application of the abovementioned tests allows for a more in-depth evaluation of aspects such as sensitivity, directionality, or comprehensiveness.
However, the results obtained in this work show that, when assessing the qualities of a method, it is not only necessary to contrast the seven parameters proposed by Hiscock and Tabrett (2010), but also to conduct the same inferential statistical tests that will then be applied to the archaeological assemblages. For example, in the specific case of Weibull distributions, the simulation of archaeological scenarios is not only ideal for understanding the limits and biases of each method, as it is done in this study, but also to help in the interpretation the reduction results themselves (e.g., Morales 2016).
On the other hand, the limitations of the SDI have led us to explore two ways of transforming the data. The log-transform of the SDI values (LogSDI) proved to be only partially effective, since the values of the Shape parameter of the Weibull distribution, although more similar to the real ones, are still higher. In addition, this transformation did not overcome the difference in scale of the values themselves, which hinders the comparisons of reduction results. Therefore, we decided to explore the rescaling min-max rescaled values assuming a 0 to 100 scale. While this scale cannot be directly equated to the percentage of volume extracted (as a normalised LogSDI value of 80 does not mean 80% reduction), this scale-independence does make direct comparison of the Weibull distributions between the different methods more feasible.
Data normalization has allowed to obtain Weibull distributions for the SDI more similar to the real data, both for the complete sample and for each of the two simulated reduction scenarios. Nevertheless, it is necessary to be cautious with this transformation for two reasons: 1) the results obtained in the subsamples by reducing the number of cores show a greater sensitivity in the normalised LogSDI data, especially regarding the Shape value, and 2) in order to make comparisons between different assemblages, it is necessary to scale the SDI data of all the assemblages under study together as if they were one. Otherwise, the possible differences between the assemblages might be camouflaged by losing resolution.
After the exploration carried out in this work, we propose the following workflow for analysing the core reduction intensity patterns in archaeological assemblages ( Figure 8): 1. Calculation of the percentage of volume extracted using the Volumetric Reconstruction Method (when possible), the PerNCO and the Scar Density Index for the cores of the assemblage. In this study, we have tested these three methods, although other methods or parameters can be used, such as the number of exploited surfaces, the number of percussion platforms, the resulting angle between the percussion platform and the flaking surface (e.g., Douglass et al. 2018;Reeves et al. 2021), or any other index. In such cases, the performance (distribution, type of correlation, etc.) of each method with respect to the real reduction data should be experimentally checked to see if any prior transformation of the data is necessary.
2. As we have discussed, in the case of SDI, it has been necessary to perform: 1) a logtransform of the SDI values and 2) a min-max scaling of the LogSDI values. The logarithmic transformation is useful for transforming asymmetric distributions with concentrations of lower values and extremely high outliers, so that the resulting curve better fits the distribution of the real data. In this way, these transformations make it possible to obtain distributions that are more similar to the real ones, as reflected in the results of the analysis of the complete assemblage and of the different sub-samples. However, it is necessary to bear in mind that this normalization must be done on all the samples to be compared. This means that, if there are two archaeological assemblages and differences among them are to be contrasted, the LogSDI data must be normalised by taking the two assemblages together, and, in the case of adding a third one later in the comparison, the data must be normalised again for the three assemblages.
3. After this transformation, the resulting distributions are evaluated by means of histograms and density curves, and then the Weibull distribution is calculated for each index. In this way, the similarities and differences of the assemblages can be quantitatively contrasted. 4. Furthermore, the analysis of correlations will allow a deeper understanding of the concordance between the different methods and, thus, a better understanding of the possible similarities and discrepancies of the results obtained. In this part of the study, it is advisable to study the relationships of the indexes with those elements that may condition the performance of the methods, such as the size of the cores and of the original blanks (Ditchfield 2016;Lombao et al. 2019). 5. Finally, the internal variability of each assemblage can be explored, analysing the differences between assemblages not only regarding their distributions. The use of multivariate analyses such as Principal Component Analysis (PCAs) allows to combine the different reduction indexes in the same analysis. In this way, the analysis of several archaeological assemblages, or different groups of interest within the same assemblage (e.g., distinct raw materials or knapping strategies, etc.) can be explored integrating multiple complementary indexes, allowing for a larger part of the archaeological record to be covered. In our case, the PCA is positively correlated with PC1 (Figure 7) since the different methods show a good correlation among them. This is especially observable through the coefficient of determination (r 2 ), whose minimum value between the different methods is 0.53. However, in any archaeological assemblage there may not be such a clear agreement between the different approaches, so it is a good tool to evaluate or confirm the agreement between different reduction methods. Furthermore, it can be a very useful element to combine reduction analyses with other technological variables, being able to integrate in the same analysis metric variables (i.e., length, width, thickness, or volume of each core), technological attributes (angles, depth of extractions, number of percussion surfaces or flaking surfaces, among many others) and reduction intensity variables. This could allow a more in-depth study of human technological behaviour and bypass the reduction stages by adopting a more quantitative perspective within the reduction continuum.
In this study we have not considered other methods, such as the one proposed by Douglass and colleagues (Douglass et al. 2018;Reeves et al. 2021), because some of the variables considered by them may be more related to the knapping strategy applied than to the reduction intensity (e.g., number of knapped surfaces or percussion platforms). However, they could be incorporated into the study depending on the characteristics of each assemblage. Similarly, this workflow is applicable to retouched tool assemblages, where it would be possible to combine methods such as 3D-ERP (Morales et al. 2015a) and GIUR (Kuhn 1990).

Conclusion
The workflow we propose in this paper achieves a fundamental objective in the study of reduction, which is the combined application of several reduction indexes to contrast the results from different methods, as suggested by Dibble (1995). In this way, by transforming the SDI data together with the application of other indexes, such as the percentage of non-cortical surface and the VRM, it is possible to study the reduction intensity patterns on cores from archaeological assemblages in a robust and statistically reliable way. Moreover, this makes it possible to complement methods that provide more detailed information but have limitations in their applicability, such as the VRM which provides information on the original volume of the blank and the reduction as a percentage of extracted volume, together with others whose information is less detailed but are easier to apply.
Finally, the use of Principal Component Analysis (PCA) to combine several indexes provide a new insight into the internal variability of assemblages, since it can offer a more precise evaluation of the reduction degree of each core. Rather than debating the pros and cons of generalized or specialized methods, future studies should focus on finding complementary proxies to supplement and cross-verify these analyses.
Además, para evaluar la sensibilidad de las curvas de distribución de Weibull, hemos submuestreado aleatoriamente el conjunto experimental a intervalos de cada 10%, lo que nos ha permitido observar la evolución de los parámetros Shape y Scale de las distribuciones Weibull tanto de los datos reales como de los distintos índices a medida que se transforma el conjunto.
Los resultados muestran que no hay grandes variaciones en los valores de Shape y Scale de los datos reales de reducción a medida que disminuye el porcentaje de núcleos en el conjunto. Sin embargo, los métodos responden de forma diferente a estos cambios en el conjunto experimental. Así el porcentaje de superficie no cortical y el Volumetric Reconstruction Method permanecen estables a lo largo de las simulaciones. Sin embargo, la normalización min-max de la transformación logarítmica del Scar Density Index muestra mayores oscilaciones en el caso del valor Shape, mientras que el valor Scale es más estable.