Alexandre Disser
email: disser.alex@gmail.com

Sylvain Bauvais

Philippe Dillmann

Transformations of the chemical signature of slag inclusions throughout experimental refining and first shaping of bloomery iron: new methodological developments 1

Keywords: Non-metallic inclusions, multivariate classification, experimental archaeometry, iron artefacts, chemical analyses, bloom smithing

 A convenient way of sorting non-metallic inclusions in iron artefacts  Experiments for assessing the chemical impact of smithing fluxes  A spatial approach of the chemical features of iron artefacts

Abstract

Archaeometallurgical experiments have been conducted on the "Plateforme d'expérimentation des arts du feu" in Melle (Deux-Sèvres, France) in order to get a better understanding of the potential alteration of the chemical signature of non-metallic inclusions in iron artefacts. This alteration may likely be caused by the use of fluxes during the smithing work. Five experimental sequences have been performed, following distinct patterns regarding the use of flux. Cross-sections made from the experimental products have sustained an as exhaustive as possible chemical characterization, resulting in the analysis of several thousands of inclusions per sample. An adapted data treatment procedure, relying on the use of multivariate analyses and geographic information systems, has been set up in order to identify the origin of each inclusion: smelting or smithing.

Introduction

Research carried out in iron archaeometallurgy has increasingly relied on the study of slag inclusions generated during the ore smelting process. Though they have been in the focus of researchers for forty years now [START_REF] Todd | Ethiopian bloomery iron and the significance of inclusion analysis in iron studies[END_REF][START_REF] Hedges | Source determination of iron currency bars through the analysis of slag inclusions[END_REF], their study has benefitted from continuous analytical developments that now allow an exhaustive acquisition of cross-sections. The petrographic and, above all, chemical characteristics of slag inclusions are induced by the material components of the smelting system, as well as the conduct of the operation [START_REF] Crew | The influence of clay and charcoal ash on bloomery slags[END_REF][START_REF] Fluzin | Apport de l'archéométrie à la restitution de la chaîne opératoire en sidérurgie ;matériaux et procédés. Etudes métallographiques, in: 6 Ème Cours d'archéologie d'Andorre ;l'obtention Du Fer Par Le Procédé Direct Du IV Ème Au XIX Ème Siècle[END_REF][START_REF] Leroy | Approche archéométrique de la proto-industrie du fer: les conditions de réduction des minerais de fer de Lorraine et la chaîne opératoire dans la sidérurgie antérieur à l'apparition du haut fourneau[END_REF][START_REF] Serneels | Archéométrie des déchets de production sidérurgique. Moyens et méthodes d'identification des différents éléments de la chaîne opératoire directe[END_REF]. These characteristics are thus used to develop a technical approach to the production of ferrous alloys (identifying the technical process, the nature of the materials used, and determining the parameters of the operation, such as the temperature or the air inflow). They are also used at a broader scale to reconstruct the socio-economic systems integrating these production activities, particularly through provenance studies (for a recent review of the literature related to this topic see [START_REF] Charlton | The last frontier in "sourcing": The hopes, constraints and future for iron provenance research[END_REF].

However, these questions cannot be addressed through the analysis of any inclusion. Depending on their mineralogy, smelting inclusions may undergo alterations during the post-smelting stages. One of these phenomena is due to a difference in the plasticity of the minerals during the heating and hammering phases of the metal, which leads to the fragmentation of some of them (Mahé-Le-Carlier and Ploquin, 1999;[START_REF] Paynter | Spinel-rich slag and slag inclusions from a bloomery smelting and smithing experiment with a sideritic ore[END_REF], thus generating new inclusions for which the original chemical features may be altered. Furthermore, the use of fluxes (sands, clays, salts...) during forging to strip or protect the metal surfaces can lead to the formation of new inclusions that are not linked to the smelting operation. It is therefore essential to ensure that the inclusions considered are characteristic of the smelting system. This implies an accurate identification of the different kinds of inclusions contained in the metal. This identification is all the more simple as the observed population is large. Previous papers have extensively discussed the issue of representativeness at the scale of the observed samples [START_REF] Charlton | Investigating the production provenance of iron artifacts with multivariate methods[END_REF][START_REF] Dillmann | Slag inclusion analyses for studying ferrous alloys employed in French medieval buildings: supply of materials and diffusion of smelting processes[END_REF].

The main purpose of this contribution is to propose some tools aimed at addressing questions commonly raised in iron archaeometallurgy, which are difficult to solve on the basis of archaeological remains, as the operating parameters are largely unknown to us. To this end, a campaign of archaeological experiments has been performed on the "Plateforme d'expérimentation des arts du feu" in Melle (Deux-Sèvres, France). It involved two professional blacksmiths, used to working in operating conditions prior to the Industrial Revolution (ventilation through bellows, charcoal as fuel, etc.). This experimentation was the opportunity to use forging fluxes of which the nature and the chemical signature were known.

Experimental archaeology for iron metallurgy has largely focused on the reconstruction of smelting systems. Forging stages, understood here as all the operations designed to obtain finished products from raw metallic masses, has only been the subject of a small number of experimental works, even though a great many questions have been raised by the multiplication of archaeometric studies on iron artefacts. In particular, the conservation of the chemical signature characteristic of the smelting system is a crucial issue. These studies have highlighted three main factors on which this conservation depends throughout the post-smelting steps: 1) the elimination of smelting-related inclusions; 2) physicochemical modifications of smelting-related inclusions; 3) the formation of new inclusions, mainly linked to the use of addings [START_REF] Blakelock | Slag inclusions in iron objects and the quest for provenance: an experiment and a case study[END_REF][START_REF] Charlton | Investigating the production provenance of iron artifacts with multivariate methods[END_REF][START_REF] Dillmann | Slag inclusion analyses for studying ferrous alloys employed in French medieval buildings: supply of materials and diffusion of smelting processes[END_REF].

In order to specify the influence of these factors on the preservation of the initial chemical signature, it is therefore essential to carry out the most exhaustive possible chemical characterisation of metallic samples taken at different stages throughout the forging process. Thus this work is based on recent developments of analytical devices commonly used in iron archaeometallurgy to develop a workflow that allows an exhaustive analysis of the chemical signature of slag inclusions trapped in ferrous products. These developments make it possible to chemically characterize almost all the inclusions visible on a cross-section, thus providing exhaustive chemical information at the sample scale. The large quantities of data generated (about fifteen variables for several thousand individuals) required the implementation of an adapted data processing protocol, relying on a customized data processing with R, and the use of GIS to facilitate the spatial display of various queries.

Experimental protocol

The raw metal used has been obtained experimentally by M. Leroy and P. Merluzzo by smelting an Aalenian oolitic ore (the "Minette Lorraine") on the "Musée de l'histoire du fer" experimental platform (Nancy-Jarville, France). This ore has been chosen because of its very particular chemical composition (high relative contents of Mg, Al, P and Ca). This peculiar chemical signature was expected to facilitate the identification of smelting-related inclusions within the experimental products. A series of five operating sequences has then been undertaken on the "Plateforme d'expérimentation des arts du feu" in Melle. It has been carried out by a professional blacksmith, who is used to forging metal obtained from the "Minette" ore, and to work with charcoal as fuel. Three crude iron blooms have been refined and shaped to bars; the two largest specimens were split into two smaller ones (1-1, 1-2, 3-1 and 3-2). All the operating sequences followed the same protocol: compacting and cleaning of the bloom; first elongation of the refining mass for bar shaping (fig. 1); bending/welding of the bar; second elongation of the bar; second bending/welding; third elongation of the bar (fig. 2). The nomenclature of the samples is as follows XPEPU15_A_B_C, where A is the number of the original bloom, B is the number of the sequence, and C is the stage at which the metal has been sampled (1: at the end of the refining stage, 2: after the first bending/welding operation, 2: after the second bending/welding operation). Experimental sequences are summarized in Table 1. In order to assess the quantities of forging inclusions introduced into the mass of metal during the operation, we divided each sequence into two phases of flux use, the first corresponding to the cleaning and compacting of the mass of raw metal, and the second corresponding to the shaping of the metal (stretching of the bar and bends/welds). During the refining stage, approximately a handful of flux has been put on the metal. While the use of fluxes during this early stage has not been documented archaeologically, we wanted to assess to what extent the porosity of the crude metal could have favoured the generation of new inclusions through pollution from the refining hearth. During the shaping of the bars and the bending/welding operations, flux has been used in such a way as to fully cover the surface of the metal. Two types of fluxes likely to have been used in a forge workshop context were used. The first is a slag collected at the Castel-Minier ironmaking workshop (Ariège, France). This slag contains fair amounts of Manganese that is intended to be used as a tracer element. The second one is a sodic glass that was obtained by crushing scraps of glass from a craft glass factory (Silycibine workshop in Arcueil, France). Both materials have been ground to a grade of about 100µm. These two materials have previously been chemically characterized by SEM-EDS. Their chemical composition is shown in fig. 3. The "Minette slag" data are given as reference for the smelting-related inclusions. The amounts given are the mean values of 11 slag samples produced by the experimental smelting of "Minette" ore. The amounts for "Castel-Minier" slag are the mean values of 17 measures made on the powder. The data for the sodic glass are the mean values of 21 measures made on the powder. Three samples were taken throughout each sequence: a sample of the raw bloom at the end of the cleaning phase; a sample of the bar before the first bending/welding; a sample of the bar after the second bending/welding. By this way, it is possible to evaluate the quantity of slag inclusions that are expelled from the metal during the process, as well as the quantity of forging inclusions that are generated.

Analytical approach

Sample preparation and micrographic exams

The samples were embedded in epoxy resin and polished (from grade 80 to a 1µm diamond suspension). Each cross-section was then subjected to a metallographic examination after acid etching (Nital with 3 vol% of HNO3 in ethanol). This first examination revealed the distribution of carburized structures, as well as features linked to the welds. A second examination of the polished section was carried out after repolishing, in order to examine the inclusions (distribution, quantity, petrography). A mosaic of pictures was acquired at each of these stages, providing comprehensive digital images of each cross-section.

Chemical characterization of the inclusions

The cross-sections were then subjected to elemental analysis by energy dispersive X-ray spectrometry coupled with scanning electron microscopy (SEM-FEG JEOL 7001-F equipped with a silicon-drift detector). The accelerating voltage was set at 15 kV and the probe current at about 8 nA. This setup allows reliable quantification for elements whose amounts are superior to a threshold comprised between 0.3 and 0.5 wt%, depending on the chemical element considered. Through this analytical device, coupling imaging and quantification of major chemical elements, several thousands to tens of thousands of entities can be characterized based on brightness values (levels of grey) obtained from backscattered-electron imaging. Metal has high values (white) while inclusions, porosities and corrosion products have low values (grey to black). Each aggregate of dark pixels (inclusion, porosity or corrosion) is analyzed and defined as an entity. Using the automated microscope motion system, each entity is spatially referenced. All image and chemical data have been primarily processed using Aztec software (Oxford Company). An automated filter has been applied on each spectrum to remove the background contribution. Quantification of chemical elements has been processed using XPP matrix correction routine (with a Phi-Rho-Z function). The following elements have been systematically quantified: O, Na, Mg, Al, Si, P, S, Cl, K, Ca, Ti, V, Cr, Mn, Fe. The results have been normalized to 100%. A spreadsheet comprising spatial, morphological and chemical information for each particle is generated at the end of this process.

Data processing procedure

The large volumes of data that are acquired have necessitated the establishment of an adapted data processing protocol, which will be described in this section. The aim is to obtain a data set that can be analyzed by means of multivariate statistics. To this end, the adopted data processing allows the removal of outliers (entities that are not inclusions or for which the analysis conditions are not good enough) and the processing of missing data (mainly related to measure below the detection limits). This procedure is transcribed in the form of a program in R, whose code is provided as Supplementary Materials. Each step described below is tagged within the script.

Dealing with detection limits

The amounts of some elements of interest (particularly Na, Mg and Mn) are below detection limits for a significant number of inclusions. Even low concentrations of these elements provide useful information for defining chemical groups of inclusions in samples. We have therefore chosen not to exclude these inclusions, as this would mean losing significant amounts of information. We therefore opted for imputation of low concentrations. We are in the case of "censored data", i.e., the low concentrations are in a range bounded by a zero concentration and the limit of detection (LoD) threshold. Several procedures are available, and have been the subject of numerous publications (a review is proposed by [START_REF] Helsel | Much Ado About Next to Nothing: Incorporating Nondetects in Science[END_REF]). We discarded the substitution of missing data with a specific concentration (e.g. half the LoD), as this causes biases in multivariate analysis (e.g. generation of linear clouds at the factorial level). We therefore opted for random generation of concentrations in the range between zero and the limit of detection. A classical approach consists in generating these contents according to a normal distribution. We rejected this strategy because it assumes that the majority of the levels would be equal to half the detection limit (modal value of a normal distribution). Therefore, we opted for a uniform distribution generation, which seems to best fit the present situation (equal probability of the grade in question being between null value and the detection limit) (Little and Rubin, 2019, pp. 67-69).

Sorting SEM-EDXS particles out

An automated procedure is then applied to eliminate aberrant entities and keep only the inclusions.

The first selection criterion is the accuracy of measurement for the most abundant inclusion compounds. Accuracy is expressed as the ratio between a measured content and its associated measurement error. The ratio used is 0.1: the measured value must therefore be at least 10 times greater than the measurement error. This ratio is calculated successively for O, Fe, Si, Al and Ca. If only one of these conditions is not met, the inclusion is eliminated from the dataset.

Taking into account different kinds of ferrous compounds

The conditions of particle analysis (image contrast, grey level thresholds, sample surface geometry, particle shape) may cause the metal matrix to contribute to a greater or lesser extent to the chemical characterization. In order to avoid this bias, a filter based on the stoichiometric balance of inclusions was used. Iron is assumed to be mainly present under valence II (forming wüstite FeO or fayalite Fe2SiO4). Each elemental concentration is converted into oxide. The sum of the concentrations expressed as oxides is therefore close to 100%, depending on whether the iron measured is of valence II. In the case of an important assay of the matrix (or metallic iron), this sum differs strongly from the theoretical sum of 100%. Different thresholds have been tested in order to avoid discarding too much valid particles and keeping biased ones. For this study, a sum of oxides lower than 80% or higher than 120% has proven satisfactory, albeit narrower intervals could be used for a stricter selection.

Defining chemical groups of inclusions

The multivariate approach to identifying chemical groups of inclusions proposed here builds on previous studies that have used major elements for reconstructing technical systems [START_REF] Charlton | Explaining the evolution of ironmaking recipes -An example from northwest Wales[END_REF] or identifying the origin of the metal [START_REF] Disser | Iron supply for the building of Metz cathedral: new methodological development for provenance studies and historical considerations[END_REF].

Raw elemental contents cannot be used, as they are highly dependent on the proportion of iron oxides within the inclusions, which is highly variable. Moreover, some elements (such as Si or Ca) are present in contents up to ten times higher than those of other elements (Na, Mg or Mn). To prevent the analyses to be over influenced by these factors, a transformation has to be applied to the data.

The various means of data transformation applied to archaeometallurgical problems have been widely discussed in several articles [START_REF] Aitchison | The statistical analysis of compositional data[END_REF][START_REF] Baxter | Log-Ratio Compositional Data Analysis In Archaeometry[END_REF][START_REF] Charlton | Investigating the production provenance of iron artifacts with multivariate methods[END_REF][START_REF] Leroy | The medieval iron market in Ariège (France). Multidisciplinary analytical approach and multivariate analyses[END_REF]. In the present study, the chemical composition is expressed in centred and reduced log-ratios. The transformed variables, referred to as Xij, were calculated using the logarithms of Si, Al and Ca contents as denominators of the ratio, as they are the most abundant elements in the inclusions (apart from O and Fe), and hence the most accurately measured elements [START_REF] Disser | Iron supply for the building of Metz cathedral: new methodological development for provenance studies and historical considerations[END_REF]. A multivariate Ascending Hierarchical Clustering treatment (hclust function from stats package) is applied to this transformed data (xij_Na, xij_Mg, xij_Al, xij_Si, xij_K, xij_Ca, xij_Mn) to define chemical groups of inclusions. The "ward" algorithm has been used to perform the clustering. Factorial analyses (Principal Component Analysis on transformed data) and matrices of biplots (plotting the elemental amounts for each couple of chemical elements) are also generated, in order to assess the relevance of the clustering.

This data processing is carried out using an R program that we have set up, the code for which is provided as Supplementary Material. The program was run in the RStudio user interface. Several graphs involved in the interpretation of the data (factorial planes, dendrograms, matrices of biplots) are generated by the script and can be recorded through the interface. The processing also provides a dataset containing all selected inclusions (comprising classification results, relative coordinates, dimensions of the inclusions and raw and transformed elemental contents).

Spatial data treatment

This dataset is integrated into a Geographic Information System managed by the QGIS software (version 3.10). The micrographic mosaic is incorporated into this system according to a custom georeferencing system, which makes it possible to examine the organization of inclusions according to their chemical characteristics. For this, we use the software's default georeferencing plugin (https://docs.qgis.org/3.10/en/docs/user_manual/plugins/core_plugins/plugins_georeferencer.html) . The mosaic is implemented by defining three to four control points. To do so, remarkable inclusions are selected within the point cloud. The transformation performed is a "thin plate spline" and the resampling method is the "nearest neighbour". The superimposition of the inclusion point clouds and the micrographic mosaic of the samples allows the groups of inclusions to be compared with the remarkable characteristics of the metal matrix, in particular the weld lines. It then becomes possible, from queries involving different kinds of information (chemical, dimensional, spatial and petrographic), to interpret each group of inclusions: forging inclusions, hercynite inclusions and smelting-related inclusions. This identification, the final stage of the analytical approach, is essential to describe the evolution of the chemical signatures of the inclusions throughout the forging process.

Results and discussions

Clustering performance

Before examining the distribution and relative amounts of the distinct inclusions classes within the samples, it was necessary to ensure that the implemented data processing procedure, from the sorting of raw data to the definition of chemical classes of inclusions by Ascending Hierarchical Clustering, made it possible to obtain a relevant classification. The criterion used for this is the distribution of the transformed variables (the Xijs) according to the class assigned at the end of the processing. A misclassification can be highlighted by the multimodal nature of the Xij distribution for a given class. We tested the unimodal character of the Xij distributions of Na, Mg, Al, Si, K, Ca and Mn for each inclusion class resulting from the clustering. The Dip test for unimodality computed by J.A.

Hartigan was used to verify this character [START_REF] Hartigan | The Dip Test of Unimodality[END_REF]. The detailed results of these tests are given as Supplementary Materials. The test was performed on 945 Xij series (7 Xij per group, 137 groups for 15 samples). The distribution was unimodal for 945 series, i.e. 98.54% of the cases. The partitioning procedure adopted was therefore effective in highlighting chemically homogeneous classes of inclusions.

The identification of the classes of inclusions will be used to describe the distribution and the evolution of the proportion of forging and smelting inclusions. The importance of the fragmentation [START_REF] Dillmann | Slag inclusion analyses for studying ferrous alloys employed in French medieval buildings: supply of materials and diffusion of smelting processes[END_REF]Mahé-Le-Carlier and Ploquin, 1999) previously highlighted in the case of the "Minette" can also be accurately assessed. Finally, the classes will be used to discuss the criteria for recognizing and interpreting forge inclusions on archaeological samples.

Interpreting inclusion classes

In order to be able to interpret each class of inclusions, their chemical features (expressed as average contents of each element of interest) were first compared with those of the reference materials: ground glass and Castel-Minier slag, and experimental smelting slag from the "Minette". The sample XPEPU15_1_1_2, obtained after using fluxes during the refining stage, the shaping of the bar and the first bending/welding, will be used as an example.

Fig. 4: Cross-section of sample 1_1_2 displaying the inclusion classes obtained through the data treatment procedure. The size of the marker is proportional to the surface of the analyzed particles.

The sampling was carried out after the first bending/welding of the bar obtained from the bloom. Since forging inclusions are generated by the reaction of the fluxes with iron oxides that form on the surface of the metal, the level of iron oxides is significantly higher than in the materials used as flux.

In order to compare the absolute contents of the chemical elements, the contents have been normalized to iron content (method described in more detail in Disser et al. 2014). Group 2 shows a high Na content and low Al and Ca contents. These chemical features are similar to those of the glass used as a forging flux in sequence 1_1, and lead to the interpretation that these inclusions are forging-related inclusions. The distribution of inclusions in the metal confirms this hypothesis, with group 2 being clearly located on the weld line (Fig. 4). Group 3 shows chemical features that are distinct from those of glass and smelting slag, and is characterized by much higher Al contents. Examination of these inclusions under metallographic microscopy shows that they consist of angular crystals of spinel identified as hercynite (Fig. 5). These crystals are commonly present in the slag resulting from the smelting of "Minette" [START_REF] Leroy | Archéologie du fer en Lorraine[END_REF]. Some of these crystals appear to be "expelled" from the smelting inclusions due to their plastic deformation while the metal is hammered. These crystals generate inclusions that will be labeled as "hercynite" inclusions.

The chemical features of groups 4, 5 and 6 are very consistent with those of the experimental "Minette" smelting slag (Fig. 6). There is little doubt that these inclusions have been generated by the smelting operation. Therefore, the inclusions from these three groups carry the smelting system's original chemical signature that can be used to determine the smelting parameters or to identify the origin of the metal.

Fig. 6: Main chemical features (mean values) of the reference materials and of the six inclusion groups defined through the data treatment procedure Inclusions forming group 1 stand out from the other groups, and do not appear related to either the forging flux or the smelting slag. These features are correlated with a very high amount of oxygen. This suggests that these inclusions contain high amounts of iron oxides, and may probably be generated partly by the oxidization of the surface of the metal while heating it. They would therefore be related to the forging process.

The approach presented through this section has been applied to each of the 15 samples taken during the forging experiments. The interpretation of the groups of inclusions will be used to support the results and discussions developed below.

Quantitative assessment of the proportions of distinct inclusion types

As said in the introduction, the formation of inclusions related to forging operations can complicate the determination of the characteristic chemical signature of the smelting system, thus hindering the identification of the source of the metal. This issue has been raised frequently, and has already been the subject of previous studies aimed at assessing the impact of forging operations on the chemical signature of ferrous alloys [START_REF] Dillmann | Slag inclusion analyses for studying ferrous alloys employed in French medieval buildings: supply of materials and diffusion of smelting processes[END_REF]. Performing extensive analyses of crosssections obtained specifically by experimental means offers the opportunity to evaluate the quantities of inclusions formed during the forging process. After interpretation of the groups of inclusions, the proportions of forging and smelting inclusions (as well as the "hercynite" ones) were calculated (Fig. 7). These proportions are expressed by the surface occupied by each interpreted type of inclusion on the polished surfaces of the samples taken at the end of each forging sequence. The These results plead for a strong conservation of the chemical signature of the smelting system in ferrous products. However, it should be kept in mind that these sequences were not intended to obtain a fine-worked metal. It can be compared to the quality of products such as hardware elements or architectural implements. Experiments aimed at obtaining finely worked products, such as blades, should be carried out to pursue these thoughts. In any case, the proportions observed (from 50 to 90 percent) for this study are large enough to suggest that even complex and fine products still contain a significant amount of smelting slag inclusions. Integrating spatial information to the protocol provides an opportunity to examine the distribution of different types of inclusions within the metal. It is of particular importance as it can develop an accurate analysis of how new inclusions are formed during post-reduction operations. A central issue in inclusion studies is to determine whether forging inclusions are exclusively related to features linked to forging (in this case, folds/welds) or are more scattered within the metal.

Although the use of fluxes during the refining phases is not archaeologically evidenced, we used this practice for sequences 1_1 and 1_2 to simulate a case where the formation of post-smelting inclusions would be maximal in a porous material. Samples 1_1_3, 1_2_3 (use of fluxes during refining stage) and 3_2_3 (no fluxes during refining stage) were compared (Fig. 8 to 10). In all three cases, the inclusions identified as related to forging did not appear randomly distributed within the metal. They constitute coherent sets, which develop parallel to the hammering planes. Most of these sets form lines that correspond to the folding/welding surfaces. The comparison of these three experimental products shows that the incorporation of fluxes during the cleaning phase (1_1_3 and

Frequency

Forging inclusion Hercynite Smelting slag inclusion 1_2_3) is negligible. In all three cases, the forging inclusions identified appear clearly associated with the use of flux during the two folding/welding stages, during the elaboration of the iron bars. The number of experimental samples is too small to be able to generalize, but these initial results suggest that the inclusions formed by forging operations are much localized within the metal, as they are overwhelmingly associated with the weld lines. Not only does this experiment provide us with crucial data on the behaviour of non-metallic inclusions in the metal during the bloom refining and forming phases, it also provides us with important information on the conduct of the processes and the behaviour of the elements as a whole. During its refining and post-smelting shaping, the metal is successively heated in an oxidizing atmosphere and a layer of iron oxide forms on the periphery of the mass. Indeed, with each pass through the fire, part of the metal is lost in the form of oxide. Through their use for removing the oxide scales from the metal surface, fluxes may contribute to prevent excessive oxidation and reduce metal loss. The experiments carried out so far efshow that between the bloom and the bar prepared to be transformed into an object, between 50 and 80% of metal is lost depending on the initial compaction of the bloom and the technical quality of the blacksmith [START_REF] Crew | The experimental production of prehistoric bar iron[END_REF][START_REF] Crew | The experimental production of bar iron[END_REF][START_REF] Dillmann | Refining of an experimental biscaïan bloom from Agorregi[END_REF][START_REF] Leblanc | Sur l'identification de la chaîne opératoire de forgeage du fer antique, associée à l'archéométrie des battitures[END_REF][START_REF] Leroy | La sidérurgie dans l'est de la Gaule. L'organisation spatiale de la production de l'âge du Fer et haut Moyen Age[END_REF] (Table 2). Within the framework of this experiment, it is possible for us to provide new data on these losses, but above all to compare the importance of these losses according to the use of fluxes, the moment of their use and their type.

% [START_REF] Crew | The experimental production of prehistoric bar iron[END_REF][START_REF] Crew | The experimental production of bar iron[END_REF][START_REF] Dillmann | Refining of an experimental biscaïan bloom from Agorregi[END_REF][START_REF] Leblanc | Sur l'identification de la chaîne opératoire de forgeage du fer antique, associée à l'archéométrie des battitures[END_REF][START_REF] Leroy | La sidérurgie dans l'est de la Gaule. L'organisation spatiale de la production de l'âge du Fer et haut Moyen Age[END_REF] Thus, during the first stage of bloom refining, the losses of 50 to 56% achieved during our experiments are quite comparable to the average losses recorded in the other reconstitution works (Table 3). However, it is notable that there is no significant loss difference whether using fluxes or not. This amount of loss seems more correlated with the size and weight of the initial blooms. For the first stage of elongation and folding, the losses, between 20 and 28%, are also very close to the average losses during the first stages of bar shaping. At this stage, a slight difference seems to appear in terms of loss depending on the use of flux or not, and on the use of glass or slag. Bloom 1 shows a loss rate of 20% using glass (sample 1_1_2) and 28% with the use of Castel-Minier slag (sample 1_2_2). Finally, the second stage of shaping and kneading shows rates ranging from 16 to 27%, which are perfectly in line with the average of previous experiments. However, using flux or not had a significant impact on the loss of metal during the last steps of the process. The rates are of 16 and 18% for the two sequences during which glass was used (respectively 1_1 and 3_2). The figures are of 24, 24 and 27% for the three sequences (1_2, 2 and 3_1) during which no flux has been used. The use of slag does not appear to have a significant protective role as it resulted in the same loss rates as when using no flux.

The use of fluxes appears to get more and more effective as the work and the cleaning of the metal progress. It is also remarkable that the choice of the type of flux plays a crucial role in protecting the metal from oxidation. Getting rid of the residual smelting slag in the bloom during the refining may require the use of more and more flux as the cleaning process progresses.

In terms of operating steps, we also tried to understand the distribution of forging inclusions in the primary cleaning steps, where the porosity of the metal is the highest and the potential trapping of fluxes is the greatest (see section 4.4). We were able to show that, in samples 1_1_3 and 1_2_3, inclusions from glass fluxes do not appear to have been generated in any significant proportion during this initial stage of work.

Thus, contrary to what had been expected, refining has turned out to be the less propitious work sequence for generating forging inclusions by using fluxes. Firstly, the porosities of the bloom do not seem to trap significant amounts of flux during this phase. Second, the use of flux itself in this phase seems questionable because of its lack of effectiveness.

Conclusion

This article is a continuation of research aimed at understanding the chemical filiation of non-metallic inclusions in objects. Since the beginning of the development of these issues, the role of flux in forge have been put forward as a limitation, or at least, a possible constraint to provenance studies based on these inclusions. It aims to provide new insights about the importance of this pollution during forging, on its recognition and on its behaviour. This contribution took profit of recent developments of analytical setups (SEM-FEG and Silicon-drift detectors) to improve the understanding of the generation and alteration of inclusions in the metal throughout the smithing process. It relies on a versatile workflow to deal with the increasing size of datasets that are generated by the chemical characterization of ancient ferrous alloys. This process produces dense point clouds whose conformation is confronted to structural features of the metal to greatly improve the interpretation of the distinct chemical inclusion groups. An experimental approach, using materials with known chemical characteristics, was used to validate the data processing method. This method can therefore be used routinely on archaeological materials for an efficient identification of slag inclusions, an essential step of provenance studies.

This study also provides further information about the efficiency of the use of fluxes throughout the refining and smithing sequences. It suggests that fluxes are of poor help for refining, and that fusible and lightly ferrous materials such as glass may be more efficient than other materials, such as the Castel-Minier smelting slag. These results will be consolidated by other archaeological experiments, focusing on the chemical filiation between the fluxes used, the formation of new inclusions and the composition of the forging slag, particularly that of the sandy-clayey slag. Compared to the archaeological data from workshop sites, this would provide a better understanding of the techniques of anti-oxidants adding and the technical and cultural choices made by blacksmiths.

Aknowledgements

The authors want to thank Thierry Tonnelier and Marion Pillone, the skilled blacksmiths who have ensured the smooth running of the experiments. The experimental session has been conducted within the frame of the Programme Collectif de Recherche "Paléométallurgie et expérimentations -Recherches sur les chaines de production des métaux aux périodes anciennes" (leader : Florian Téreygeol) conducted at the experimental platform at Melle (France). They also want to thank all the people who participated in this experimental session and in particular Marc Leroy, Paul Merluzzo, Enrique Véga and Jean-Charles Méaudre.

Fig. 1 :

 1 Fig. 1: Sequences 1_1_1 and 1_2_1 of bloom refining and shaping of the two bars (hot cutting of the refining bloom into two bars.). Reading from left to right and top to bottom.

Fig. 2 :

 2 Fig. 2: Sequences 1_1_2 and 1_1_3 of elongating the bar, first folding and welding, elongating the bar and second folding and welding. Reading from left to right and top to bottom.

Fig. 3 :

 3 Fig. 3: Main chemical features of the experimental smelting slag (Minette) and the two fluxes (Castel-Minier slag and glass)

Fig. 5 :

 5 Fig. 5: Top: Secondary electron image showing an inclusion with hercynite crystals; bottom: EDS spectrum acquired on the highlighted crystal

  smelting inclusions are dominant in the final sample of all sequences, regardless how frequently flux addings have been used. The lowest proportion remains slightly above 50% (sequence 1_1).

Fig. 7 :

 7 Fig. 7: Frequency of the three types of inclusion in the final product of each experimental sequence

Fig. 8 :

 8 Fig. 8: Distribution of the three inclusion types on the cross section of sample 1_1_3

Table 2 :

 2 Summary of metal loss proportions by work phase based on

			of losses
		min.	max.	av.
	Bloom refining	39	74	52
	2nd refining	19	44	28
	Bar elaboration	4	45	19
	Bloom and 2nd refining	45	80	63
	2nd refining and bar elaboration	41	64	52
	Complete refining	52	83	73
	Object elaboration	4	27	13

Table 3 :

 3 Mass of original blooms, intermediary and final metal pieces, with corresponding loss rates

		Bloom 1 (glass)	Bloom 2 (glass)	Bloom 3 (none)
	Initial mass in g	2222	1817	4284
	After Step 1 in g (bloom refining)	1018	917	1887
	Loss		54%	50%		56%
		Mass 1-1 (glass)	Mass 1-2 (CM slag)	mass 2 (none)	mass 3-1 (none)	mass 3-2 (glass)
	Mass after 1 st sampling in g	484	534	841	863	951
	After Step 2 in g (1 st folding and
	welding)	385	386	652	639	756
	Loss	20%	28%	22%	26%	21%
	Mass after 2 nd sampling in g	341	345	608	609	683
	After Step 3 in g (2 nd folding and
	welding)	286	262	460	447	560
	Loss	16%	24%	24%	27%	18%

This paper is part of a Special Issue on the Contributions of Experimental Archaeology to Excavation and Material Studies, guest edited by Frederik W. Rademakers, Georges Verly, Florian Téreygeol and Johannes Auenmüller.


Fred Gonsalves
email: fred-michael.gonsalves@imt-atlantique.fr

Bastien Pasdeloup

Romain Billot

Patrick Meyer

Arnaud Jacques

Matthieu Lorang

New Insights into the Propulsion Power Prediction of Cruise Ships

Keywords: propulsion power prediction, methodology, neural networks

Ship propulsion is the largest consumer of energyand by extension fuel -on cruise ships. Improving its efficiency is thus an important aspect of energy management, both for environmental and economic reasons. Various approaches have been detailed in the literature for improving propulsion efficiency, ranging from optimal voyage planning to prediction of propulsion power or fuel consumption using Machine Learning algorithms, trained on high frequency sensor data. On this latter topic, the approaches typically involve a series of data transformations and time-aggregations (windowing), followed by shuffling and separation of data points into train and validation sets. However, this approach leads to very similar data in the train and validation sets, preventing trained models to generalize well on future ship voyages. In this article we highlight methodological issues and give insights on how to tackle them to train models that focus on optimizing generalizability, especially predictive accuracy on unseen future test sets. We present a temporal approach to splitting data into train, validation and test sets. We perform our analysis using simple multilayer perceptron architectures, of distinct dimensions. Our study concludes that smaller/simpler models, trained on temporal-split data have a lower error when predicting on unseen future test data, compared to larger models and usage of shuffle-split datasets, while also providing better confidence in model accuracy, due to reduced discrepancy between obtained validation and test errors.

I. PROBLEM STATEMENT

A. Context and Challenges

Global shipping accounts for roughly 2.8% of global carbon emissions, according to the International Maritime Organization 2014 Green House Gas study [START_REF] Organization | Third IMO GHG Study 2014[END_REF]. The current objective of the organization is to reduce shipping carbon intensity by 40% by 2030 and 50% total (70% intensity) by 2050 compared to the 2008 baseline [START_REF] Organization | Third IMO GHG Study 2014[END_REF]. It specifies evaluation indexes, as well as milestones to be achieved, based on the type and size of the ship. Emissions of other pollutant chemicals such as sulfur are also increasingly being tightened by restrictions internationally and in protected zones [START_REF]IMO2020 -cutting sulphur oxide emissions[END_REF].

Shipping companies and shipyards alike are thus under pressure to reduce fuel and energy consumption and, consequently, carbon emissions to comply with these relatively recent objectives and evolving regulations. In this context, cruise lines play a crucial role in reducing their greenhouse gas and chemical emissions owing to size and number of passengers. Cruise lines and shipyards work in tandem to innovate fuel-saving solutions involving all phases from design to operation of the ship. While a cruise ship can be compared to a small city in the number and variety of energy consumers on board (e.g., air conditioning, hotel, water production, lighting, entertainment infrastructure), the major fuel consumer of a cruise ship is the propulsion system.

The propulsion of a cruise ship accounts for roughly 40-60% of the ship's total energy consumption [START_REF] Baldi | Energy and exergy analysis of a cruise ship[END_REF]. The performance of the ship in terms of its energy efficiency is thus highly dependent on this system, and efforts are being made in various directions to improve it. Energy efficiency improvement strategies can roughly be classified as either design or operational. Design improvements may include hull form, propeller design, energy conversion efficiency, use of greener fuels and alternative energy sources (e.g., liquefied natural gas, methanol), or supplementary equipment (e.g., air lubrication system). From an operational perspective, measures such as speed and trim optimization, best route prediction, or regular hull maintenance can lead to important energy savings. With respect to the latter point, hull fouling -i.e., accumulation of marine growth on the hull -could indeed lead to an increasing performance deterioration and hence an increase in energy consumption [START_REF] Laurie | Machine learning for shaft power prediction and analysis of fouling related performance deterioration[END_REF].

Analysis and evaluation of energy-saving measures has been boosted by the improvement of data collection, storage and processing technologies. The literature today contains various examples of monitoring tools capable of collecting high frequency data from various sensors and analyses are being carried out over several months of data. Among these tools, Machine Learning (ML) has gained a lot of interest, especially models based on Artificial Neural Networks (ANNs). However, it appears that most works from the field of propulsion power prediction suffer from methodological issues, which is the main message of this article.

B. Related Work

In recent years, various studies have analyzed the use of ML algorithms applied to ship propulsion. The problem is typically stated as a prediction problem with the target variables being propulsion power [START_REF] Pedersen | Prediction of full-scale propulsion power using artificial neural networks[END_REF], [START_REF] Petersen | Statistical modelling for ship propulsion efficiency[END_REF], [START_REF] Theodoropoulos | Evaluation of different deep-learning models for the prediction of a ship's propulsion power[END_REF], fuel consumption [START_REF] Kim | Development of a fuel consumption prediction model based on machine learning using ship in-service data[END_REF], [START_REF] Uyanık | Machine learning approach to ship fuel consumption: A case of container vessel[END_REF], or speed [START_REF] Mao | Statistical models for the speed prediction of a container ship[END_REF], [START_REF] Brandsaeter | Ship speed prediction based on full scale sensor measurements of shaft thrust and environmental conditions[END_REF] and the input variables being operational parameters such as ship trim, draft, speed, and external conditions such as wind and sea state. The datasets used are either noon reports (data sheets prepared on a daily basis by a ship's chief engineer, consisting of average values of the abovementioned parameters, as well as the total fuel consumption and other information) or ship sensor data, sampled at varying frequencies, weather forecast or hindcast data. In the field of propulsion power prediction, the first ML-based approachbased on sensor data -was introduced in 2009 [START_REF] Pedersen | Prediction of full-scale propulsion power using artificial neural networks[END_REF]. However, the use of such data has only seen an increase starting around 2017, and the number of papers on the topic yet remains limited.

Although data can be viewed as signals or time-series, a common simplification of the problem -assumed by the majority of articles to date -consists in considering data as a collection of independent observations. While this certainly reduces the accuracy of predictions, it allows to evaluate first the predictive powers of features without considering their dynamics. In this article, we choose to make the same simplification, as we want to focus on methodological issues in the literature. Considering time series, though, is a clear line for future work.

The standard approach taken in the literature for preprocessing data before training and testing ML models thus consists of the following steps:

1) applying certain filters [START_REF] Themelis | Data acquisition and processing techniques for a novel performance monitoring system based on KPIs[END_REF]; 2) extracting features;

3) performing appropriate aggregations (e.g., mean, variance, derivative) of the input and target variables. The data points are then split into train, validation, and/or test sets. Although, performing this split is approached in various ways in the literature. The most common one is the shuffle-split approach, in which the processed data points are randomly shuffled and then allocated to these subsets.

Pedersen et al. [START_REF] Pedersen | Prediction of full-scale propulsion power using artificial neural networks[END_REF] and Petersen et al. [START_REF] Petersen | Statistical modelling for ship propulsion efficiency[END_REF] are the first studies focusing on ANNs and sensor data, collected at a high frequency (up to 1Hz in [START_REF] Pedersen | Prediction of full-scale propulsion power using artificial neural networks[END_REF]). Data is then aggregated using a fixed window length. The architecture implemented in both cases is a MultiLayer Perceptron (MLP), with hyperparameters selected using a parameter sweep and k-Fold Cross-Validation (k-Fold CV). Petersen et al. [START_REF] Petersen | Statistical modelling for ship propulsion efficiency[END_REF] propose two methods for shuffling data: 1) shuffling individual 10-minute windows; 2) shuffling trips. The authors report that shuffling trips leads to a slightly higher prediction error, but that this method avoids what they refer to as cross-talk, whereby the data in the train and test sets are likely to be very similar. This is very reasonable, since the dynamics of large moving ships are long-term and consecutive ten-minute aggregated windows of data are not likely to be very different from each other.

Drawing on the conclusions of [START_REF] Pedersen | Prediction of full-scale propulsion power using artificial neural networks[END_REF], Du et al. [START_REF] Du | Two-phase optimal solutions for ship speed and trim optimization over a voyage using voyage report data[END_REF] use a MLP architecture on noon report data. The goal of their approach is to quantify the impact of sailing speed, trim, displacement, weather and sea conditions on a ship's fuel consumption rate and then to optimize speed and trim for each segment of a voyage in order to minimize the total fuel consumption.

Kim et al. [START_REF] Kim | Development of a fuel consumption prediction model based on machine learning using ship in-service data[END_REF] study the use of MLP and multiple linear regression models to predict the fuel oil consumption of a container ship using 6 months of sensor data. Training of the MLP model is performed using train/validation sets generated by a shuffle-split approach. The authors then test the model on an independent voyage not used at all during training. Similar settings are also used in the work of Farag et. al [START_REF] Farag | The development of a ship performance model in varying operating conditions based on ann and regression techniques[END_REF], in which they train an ANN on a first voyage and report test errors on a distinct second one, for an oil tanker.

Gkerekos et al. [START_REF] Gkerekos | Machine learning models for predicting ship main engine fuel oil consumption: A comparative study[END_REF] evaluate the use of various ML algorithms to predict fuel oil consumption using both noon report data and sensor data collected using an automatic data logging and monitoring system, sampled every hour. The algorithms include, among others, linear regression, tree-and forest-based models, k-Nearest Neighbors (k-NN) and MLPs. Separate models are trained for each type of dataset. The noon report dataset contains two and a half years of data, while the other contains only three months of data.Each dataset is then split into train and test subsets (80%-20%) and k-Fold CV is performed on the train set.

Uyanik et al. [START_REF] Uyanık | Machine learning approach to ship fuel consumption: A case of container vessel[END_REF] also perform a comparison of various ML models for the prediction of fuel oil consumption. Here, data comes from noon reports and the target variable is the total fuel oil consumption over the voyage. The methodology used is similar to the shuffle-split approach. Though, no aggregations needs to be performed as the information is already aggregated in the noon reports. Grid-search and k-Fold CV are used to find optimal hyperparameters for the different ML models.

Bes ¸ikc ¸i et al. [START_REF] Bes ¸ikc ¸i | An artificial neural network based decision support system for energy efficient ship operations[END_REF] implement a MLP architecture to predict fuel oil consumption from noon report and weather data. The data points are shuffled and split into train and validation sets (70%-30%). Hyperparameters for the MLP were selected using trial-and-error.

Theodoropoulos et al. [START_REF] Theodoropoulos | Evaluation of different deep-learning models for the prediction of a ship's propulsion power[END_REF] evaluate the use of two ANNs for the prediction of ship propulsion power using sensor data. First, a MLP architecture similar to the previous literature is used. Data is preprocessed using various filtering and outlier removal techniques. Data is then smoothed using a simple moving average over 5-minute windows. Training is performed following the shuffle-split approach and using k-Fold CV. The final model is then used in a simulated operational scenario where the first 80% of the data is used for training and the final 20% data for testing. The article evaluates a second neural network model for propulsion prediction using Long Short-Term Memory (LSTM) architecture. However, as mentioned earlier in the section, we will focus on models not incorporating the temporal aspect of the data.

Finally, Laurie et al. [START_REF] Laurie | Machine learning for shaft power prediction and analysis of fouling related performance deterioration[END_REF] tackle the subject of propulsion power prediction using ML models. Their study evaluates various ML models, trained with sensor data (10 second sampling rate) and weather data. The authors provide a model to analyze the effect of hull fouling on the propulsion performance of the ship. They use a shuffle-split approach after aggregating over 10 minute intervals, followed by hyperparameter selection using k-Fold CV. Laurie et al. point out that a k-NN algorithm achieves the second best performance with a Mean Absolute Percentage Error (MAPE) only 0.07% lower than the best model -a random forest.

C. Knowledge gaps and contributions of the paper

The state of the art raises several issues regarding a potential implementation of prediction models in operational conditions:

1) The cross-talk problem [START_REF] Petersen | Statistical modelling for ship propulsion efficiency[END_REF]: Separating high frequency, sensor-acquired data into train and validation subsets after shuffling would always lead to very close analogs in the train set for every validation data point (see Figure 1). An issue thus arises when using this shuffle-split approach for validating the predictive accuracy of the model: although the validation error reported for a model may be very low, it is possible that the model is strongly overfitting to training data, leading to high accuracy errors in testing conditions.

time

Step 1: shuffle

Shuffle

Cross-validation Test

Step 2: split + cross-validation + test

Fig. 1. The cross-talk problem

An indication of the existence of this cross-talk issue is the performance of the k-NN algorithm in evaluation papers [START_REF] Gkerekos | Machine learning models for predicting ship main engine fuel oil consumption: A comparative study[END_REF], [START_REF] Laurie | Machine learning for shaft power prediction and analysis of fouling related performance deterioration[END_REF]. In particular, Laurie et al. point out that for any data point in the validation/test sets, there exist very close analogs in the train set, when data is sampled at a high frequency [START_REF] Laurie | Machine learning for shaft power prediction and analysis of fouling related performance deterioration[END_REF]. Figure 2 illustrates this phenomenon. 2) Overlapping windows: The problem of cross-talk is further exacerbated when overlapping windows are used for aggregation, such as in the case of simple moving averages [START_REF] Theodoropoulos | Evaluation of different deep-learning models for the prediction of a ship's propulsion power[END_REF]. Such an operation, used in conjunction with a shuffle-split approach, causes train and validation datasets to share information, leading to trained models not being able to learn generalizable characteristics of the data. Thus, validation accuracy reported may be misleading if the goal is to implement models for predicting ship behavior for future voyages. 3) Parsimony: the literature review exhibits already a wide range of ML models, from simple (e.g., multiple linear regression) to complex (e.g., recurrent neural networks) approaches. Laurie et al. [START_REF] Laurie | Machine learning for shaft power prediction and analysis of fouling related performance deterioration[END_REF] suggest that parsimonious models such as k-NN can outperform most refined techniques for the prediction task. Therefore, one gap is related to a rigorous methodology to compare models' complexity, for instance by comparing small vs. large architectures within the same models family. As an example, a small MLP performing at least as well as a large model on future unseen data would indicate overfitting in the selection of the larger one. Moreover, overfitting could be detected when a model performs much worse on future, unseen test data than could be reasonably expected. Concretely, the test error is much higher than the validation error.

The goal of this paper is NOT to propose a new comparison of ML methods for propulsion power prediction. The purpose is rather to come up with a simple methodological contribution that addresses the above-mentioned knowledge gaps. The research questions we would like to tackle are three-fold: Q1) Given a propulsion power prediction task, what is the most appropriate splitting approach (train/validation/test) to obtain a model leading to similar validation and test errors (i.e., providing more confidence for implementation in operational conditions)? Q2) Given a fixed size of a model, which splitting approach gives the minimum test error (i.e., which approach generalizes better)? Q3) Given a specific splitting approach, does a small model perform as well on a future test as a large model (i.e., principle of parsimony)?

With respect to these research questions, the main contributions of this paper can be summarized as follows:

• Methodological framework: We propose to compare two methodologies to train and validate a machine learning model for prediction of propulsion power using 27 months of high-frequency sensor data. The first of these is the shuffle-split approach, whereby data for both train and validation datasets comes from exactly the same timeframe (e.g., 1 year). In the second approach, data is split temporally (we refer to this as the temporal-split approach); the first part of the dataset is used for training, the second part for validation. • Models dimension analysis: Based on a MLP architecture, we also analyze the dimensions of the models obtained using the two approaches, and propose a crosscomparison of small vs. large MLP configurations in both splitting strategies (shuffle vs. temporal split). • Realistic evaluation: The best models obtained from the two approaches are then tested on an independent future test set, unseen by either model during the training process. This allows us to clearly see the gaps between validation and test errors (hence assessing overfitting) and also highlight the degradation of the performance when the prediction horizon increases.

The paper is organized as follows. In Section II, we present the global methodological framework and the splitting strategies. Section III is related to the data presentation, including preprocessing, feature engineering and windowing. Results are presented in section IV. Section V concludes the paper with a discussion.

II. SYSTEMATIC METHODOLOGY

To answer the three research questions Q1, Q2, Q3 presented in Section I-C, we propose a global methodological framework that enables a comparison of the two splitting strategies for training and validation, as well as a crosscomparison of small vs. large MLP architectures. The process starts with the selection of a machine learning framework.

A. Choice of the ML framework

A MLP architecture was selected as the ML framework to be used in this study. This choice was made for several reasons:

1) This approach is majoritary in the literature on propulsion power/fuel consumption prediction; 2) The dimensions of two MLP models can easily be compared in terms of number of weights to train, related to number of layers and layer size; 3) Existing work tends to use large MLPs (e.g., [START_REF] Theodoropoulos | Evaluation of different deep-learning models for the prediction of a ship's propulsion power[END_REF]). We want to explore the use of smaller networks and their impact on prediction. MLP can have various architectures, defined by hyperparameters such as the number of layers and nodes (neurons), or choices of activation functions. While there is a contemporary dominance of neural networks and deep learning in the ML community, the goal of this paper is not to come up with the best architecture for propulsion power estimation -e.g., including convolutional layers and temporality (e.g., LSTMs). Hence, we proceed to simple MLP architectures and rather maintain focus on the learning process. For simplicity, we use an architecture where the number of nodes per layer is constant, and find optimal hyperparameters using a grid search procedure.

B. Global methodology

Given an appropriate ML framework, Figure 3 presents the flowchart diagram of the different steps performed during experimentation.

As shown in Figure 3, the systematic methodology is composed of the following blocks:

Hyperparameters

Grid (small

) (B) (C) (D) (E) (F) (G) (H) (I)
Fig. 3. Global methodological framework (H, I) Cross-training: the purpose is to use the optimal hyperparameters obtained using one approach (e.g., shufflesplit) to train a model using the dataset of the other approach (e.g., temporal-split), yielding to two extra models (described in Section II-D).

All four models are evaluated on an independent future test set, unseen by the models during the training process. 1) Shuffle-split approach: In this first approach, all raw training data is used both for training and validation. This is performed by shuffling the data points after aggregation, as described in Figure 4; 2) Temporal-split approach: In the second approach, data used for training and validation are selected using a temporal split as described in Figure 4. Train and validation set sizes are set the same for both approaches. The temporal split was applied assigning the first 12 months of data for training and the next 3 months for validation.

C. Shuffle vs. Temporal Split

D. Models and evaluation

We now describe the models that we obtain using the two splitting approaches -labeled (F) and (G) in Figure 3 -as well as the models obtained obtained after cross-training -labeled (H) and (I) in Figure 3 -. In terms of the MLP learning process, we use settings described in Table I Number of layers and number of nodes per layer were defined experimentally through grid search to find optimal hyperparameter values. We have used two distinct grids of possible hyperparameter values, deemed large (described in Table II, used in step (F) in Figure 3) and small (described in Table III, used in step (G) in Figure 3). The use of these two grids aims at addressing research question Q3 on the principle of parsimony. Values of the large grid have been chosen to be in the same range as those used in [START_REF] Theodoropoulos | Evaluation of different deep-learning models for the prediction of a ship's propulsion power[END_REF].

Parameter

Possible values Layers 10, 15, 20 Nodes per layer 50, 100, 200, 300, 400

TABLE II HYPERPARAMETERS GRID FOR THE LARGE MODEL

For each combination of layers and nodes, k-Fold CV was used to train 4 models (k = 4) to account for model initialization and batches variability. The performance of the parameters is evaluated using the average validation loss over the k models obtained per parameter. Models resulting from the training procedures using grid searches are labeled M SL and M T S in Figure 3 1 . Comparison between these two models will give us an idea of whether or not overfitting is occurring.

Next, in order to test the effect of model size independently from the splitting approach, we interchange the optimal hyperparameters found during model building and the train/validation datasets used to obtain them (step referred to as cross-training in Section II-B. This produces two more models, M SS and M T L , as shown in Figure 3.

The evaluation metric used is the Mean Absolute Percentage Error (MAPE) [START_REF] Makridakis | Accuracy measures: theoretical and practical concerns[END_REF]. The network takes as input an array of sensor readings and weather forecast data and produces as output an estimated propulsion power. The network is then trained using the MAPE as the loss function to minimize.

III. DATA PRESENTATION

In this section we discuss the data that was used for the experiments, as well as the various preprocessing steps that were taken to prepare the prediction task. For confidentiality reasons, information allowing a possible identification of the ship (e.g., absolute speed and propulsion power, geographical coordinates) is not described.

A. Data description

The dataset comes from a large cruise ship in operation. Available data consist of a few dozen variables, however, only 20 of these variables with a known effect on the overall ship resistance were selected. Two primary types of data are present in the set:

1) Sensor data: Data coming from sensors and systems installed onboard the ship. This data is used in the automation and control systems of the ship. Data is first stored on board and then transferred to the data provider.

Sensor data selected includes speed (over ground and through water), drafts (fore and aft), acceleration (over ground and through water), number of stabilizers used, air and sea temperatures, thruster power, and electrical propulsion power. 2) Weather forecast data: This data is acquired from a third party meteorological data provider. All interpolations are performed by the meteorological data provider, and information regarding how these interpolations are performed is not, at this time, known. Information collected includes relative wind (longitudinal and transversal speed) and wave (period, angle and height).

The dataset collected contains 27 months of data. The sensor data is sampled every 30 seconds. Figure 5 shows a sample time series of the target variable, propulsion power. Accordingly, the weather data also contains data points at a 30-second sampling rate -however, as mentioned earlier, this data is interpolated spatially and temporally by the meteorological data provider as a function of the ship's geographical coordinates and time. Weather forecast data was used in this study as the model aims at being used for propulsion power forecasting purposes. Air and sea temperatures were taken from the sensor readings as they were found to be very highly correlated with the weather forecast data.

The data collected contains temporal gaps and noise issued from possibly varied sources (e.g., sensor error, extrapolation error). This requires a few basic filters and preprocessing steps to be applied before the dataset can be used.

B. Preprocessing and feature engineering

As the goal of this study is not to find the optimal preprocessing steps or features that provide the best model, but rather to focus on the implications of the methodology used, only simple operations were performed for cleaning and preparing the data:

• Feature transformations: Two types of feature transformations were used. For variables that are supposed to contain integer values, a rounding was applied if the values were reported with decimal points. For example the number of stabilizers used should contain values 0, 1, or 2. For variables containing angles represented in degrees (0 -360), the sine and cosine of the angles (in radians) were extracted and used as input features; • Filters: Two types of filters were used in the processing of the 30-second raw data: 1) Operational filters: data was filtered to remove outliers concerning the speed and propulsion power of the ship. This includes setting a minimum and maximum value for the variables concerned. Data points having thruster power greater than 0 were also eliminated from the set, as thrusters are only used during maneuvers and not at sea. These filters were applied to keep only data representative of the ship's time spent at sea and to remove periods spent at berth, maneuvers and outliers; 2) Noise filters: after visualizing the various data, it was found that certain sensors were giving erratic values at certain times. These data points were manually filtered out.

C. Windowing and aggregation

Windowing and aggregation has been performed in all of the studies dealing with high frequency sensor data. The reasons for these are to remove potential noise from the sensors, and also to extract features relevant to the dynamics of the system. Studies have applied various window sizes ranging from 3 to 15 minutes depending on the dataset and application. For this study, we have used a window size of 5 minutes.

Windowing produces an M ×N matrix of the input variables (where M is the number of samples per window and N is the number of input variables) and an N × 1 vector of propulsion power. The window is rejected if any value -either in the input matrix or in the output vector -is missing. This is done in order to not deal with data imputation and its effects.

If the window has not been dropped, the data in the window is aggregated by taking the average over the time dimension, producing an input vector of dimension 1 × N and an output vector of dimension 1 × 1.

IV. RESULTS

Table IV presents the final hyperparameters selected for the two splitting approaches presented in Section II-C, i.e., the output of (F) and (G) in Figure 3. These hyperparameters represent the models showing the smallest validation loss in their respective approaches.

Using the shuffle-split approach and the large grid, optimal hyperparameters found are the largest possible values from Table II (20 layers and 400 nodes per layer). The temporalsplit method, trained with the small grid, produces a smaller network, with optimal values being intermediate in tested values from Table III According to the proposed methodology, the hyperparameters of M SL and M T S have been used on the datasets built by the two other splitting approaches, temporal and shuffle, respectively, leading to two extra models M T L and M SS . Table V With respect to the research question Q1, we notice that the two models applied to a shuffle-split data set exhibit the lowest validation losses, respectively 0.75 et 2.75 for M SL and M SS , and yet are associated with the highest test losses (6.33 and 5.52). This result demonstrates the presence of the cross-talk problem between the train and validation sets, and the degradation of the generalization performances when the test set is appropriate. On the other hand, the temporal split allows getting more confident validation losses, compared to test losses: 4.45 vs. 5.41 for M T L , 4.35 vs. 5.14 for M T S .

The aim of research question Q2 is to compare which splitting approach gives us the minimum test error, given a fixed model size. Figure 6 reports the error distributions of the four models for the 12 months of test data. Comparing large size models, Figure 6 shows that M T L outperforms M SL with lower percentiles and less variability. The same conclusion lies for smaller architectures, as we can see the dominance of M T S over M SS with a comparable interquartile range. Thus, temporal-split enables a better generalization with unseen data. Regarding the research question Q3, whose goal was to investigate the models' dimensions, Figure 6 indicates that, for a given splitting approach, a small architecture is in nearly all cases slightly better than a large one. This can be seen by comparing M SL vs. M SS , as well as M T L vs. M T S . The significance of the differences between pair-to-pair distributions have been assessed through hypothesis testing (Fisher-Snedecor procedure [START_REF] Johnston | Econometric methods[END_REF]) with a 99% confidence level.

A last category of results deals with the deterioration of the performances as the prediction horizon increases. Figure 7 shows the evolution of the monthly-averaged MAPE for the two optimal model of each splitting strategy, namely M SL and M T S . Test month index The best model obtained from the temporal-split, M T S presents a later and lower performance decrease as the prediction horizon decreases. In addition to a larger degradation from the first months, M SL exhibits a larger confidence interval, hence more variance. This trend is confirmed by Figure 8 which displays monthly APE distributions instead of single monthly averages. These last two figures also show that after a given amount of time, the two models tend to harmonize towards higher error rates, which could be explained by the presence of additional, time-evolving factors such as hull fouling.

V. CONCLUSIONS AND PERSPECTIVES

The aim of this paper was to provide new methodological insights into the problem of propulsion power prediction. After the identification of several knowledge gaps from the state of the art, we have proposed a systematic methodology in order to assess appropriate splitting strategies and model complexities for the task of predicting the propulsion power of cruise ships on unseen future data. The main conclusions can be summarized as follows:

• A shuffle-split including the test set, as shown in Figure 1, leads to the cross-talk problem, a methodological flaw yielding an overestimation of the models performances (see Figure 2); • Given an appropriate independent test set of unseen future data, a training methodology based on a temporal-split training and validation historical data set always outperforms a shuffle-split cross-validation. A shuffle split underestimates the validation error by overfitting during the training phase; • For a given splitting approach, within the MLP framework, parsimonious models are equal or better than most refined architectures, with less variance in the test errors.

Several issues merit discussion:

1) Explanation vs. Prediction: A shuffle-split methodology on high frequency sensor data would be appropriate in situations where future prediction is not the final goal. If the goal is to analyze the effects of a certain variable on a historical dataset (reanalysis), such an approach may be appropriate. In these cases, ML is employed as an explanatory tool rather than a prediction tool, and could possibly be applied to evaluating factors affecting the propulsion efficiency (e.g., [START_REF] Laurie | Machine learning for shaft power prediction and analysis of fouling related performance deterioration[END_REF]) or evaluating the efficacy of certain equipment; 2) Even though the error increases for data outside the window of the training data, the temporal split model still performs with lower error at a 6-months horizon. This suggests that, even though we might expect the model to perform worse with time as new data is encountered and hull fouling may occur, the model trained using temporal split has generalized better than its counter part; 3) Sensitivity analysis on the effect of window's size on overfitting: It is likely that larger windows (e.g., 30, 60, 90 minutes) would be able to overcome cross-talk. An extreme case would be to shuffle entire trip segments (port-to-port) when enough data is available.

Perspective work to improve on the understanding and prediction of propulsion power for cruise ships could include a combination of supervised and unsupervised algorithms, and the inclusion of physical knowledge to the problem. Within the context of supervised learning, data could be segmented by trips and the trips could be shuffled and split into train, validation, and test sets. This would likely enable the learning of a more generalized model and avoid the cross-talk problem. Furthermore, temporal aspect can be considered a) using timeseries based approaches (e.g., ARIMA models, LSTM); and b) by including time passed since hull/propeller cleaning as a feature to encapsulate the biofouling effect.

Unsupervised methods could be used to classify data based on power consumption (similar to [START_REF] Brodić | Detecting of the extremely low frequency magnetic field ranges for laptop in normal operating condition or under stress[END_REF], [START_REF] Ruiz | A time-series clustering methodology for knowledge extraction in energy consumption data[END_REF]) or power consumption deviation from an expected (theoretical) value. The latter could provide insights on factors affecting deviation from -for example -a physical model, and could be used to build a hybrid predictive model. Moreover, if entire trips are segmented, clustering could be performed at the trip-level considering consumed energy instead of power.

Fig. 2 .

 2 Fig. 2. Shuffled data

  (A) Data preprocessing: consisting of feature transformations and filters presented in Section III-B; (C, D) Data splitting: divide the dataset into train, validation and test sets according to the two possible approaches: shuffle-split and temporal-split, described in Section II-C; (B, E) Windowing and aggregation: described in Section III-C.Depending on the selection strategy (shuffle vs. temporal), this step comes before or after the split; (F, G) Train model (grid search): for both temporal-and shuffle-splits, the optimal models and their hyperparameters are found;

Figure 4

 4 Figure 4 depicts the two strategies taken for model training and evaluation -shuffle-split and temporal-split -for preparing the train and validation sets.

Fig. 4 .

 4 Fig. 4. Comparison of shuffle-split and temporal-split approaches for train and validation sets generation. Note that the test set remains identical.

Fig. 5 .

 5 Fig. 5. Sample propulsion power time series (values are normalized for anonymity)

Fig. 6 .

 6 Fig. 6. Boxplot distributions of the Absolute Percentage Error for the 4 models

Fig. 7 .

 7 Fig. 7. Monthly Test MAPE for models M SL and M T S

Fig. 8 .

 8 Fig. 8. Monthly Test APE for models M SB and M T S

  .

	Parameter	Value
	Batch Size	128
	Optimizer	Adam [17]
	Learning Rate	5e-4
	Max Epochs	50
	Loss	Mean Absolute Percentage Error (MAPE)
		TABLE I
		SHARED HYPERPARAMETERS

  (3 layers and 40 nodes per layer).

	Model Layers	Nodes per layer
	M SL	20	400
	M T S	3	30
		TABLE IV
	GRID SEARCH RESULTS: OPTIMAL HYPERPARAMETERS FOUND DURING
	STEPS (F) AND (G) IN FIGURE 3

  reports the validation and test losses of the four models.

	Model Val. Loss	Test Loss	Difference
	M SL	0.75	6.33	88.3%
	M T L	4.45	5.41	17.7%
	M SS	2.75	5.52	50%
	M T S	4.35	5.14	15.2%
		TABLE V
	COMPARISON OF VALIDATION AND TEST LOSSES FOR THE FOUR MODELS

Models are identified as Mxy, where x ∈ {S, T } denotes the use of the Shuffleor Temporal-split dataset, and y ∈ {S, L} denotes the use of hyperparameters coming from the Small or Large grid for training.