A Machine Learning Model for Prediction of Marine Icing

Icing on ships and offshore structures occurs from sea spray icing sourced from seawater and atmospheric icing due to precipitation in the form of snow, rain, fog, etc. Sea spray is either wind-generated, where seawater droplets from breaking waves are carried to the structure due to wind, or impact generated, where the sea spray is a result of the impact of the ship or structure with waves. Sea spray icing accounts for 80–90% of all offshore icing incidents [1]. About 40 fishermen are reported to have died in icing-related incidents in the 1960s [2]. Impact-generated sea spray is regarded as the most important water source in dangerous icing events [3–7]. Icing also affects offshore wind turbines, but is limited to lower levels of the structure and blade tips during azimuthal downward position [8]. Similarly, for other structures like oil rigs, sea spray icing is dominant in the lower parts where sea spray droplets are easily impacted, whereas icing in the upper regions is usually atmospheric icing [9]. Models for predicting sea spray icing are crucial for the optimal design of ships and structures operating in cold climate regions [10].

The 1960s–1980s saw extensive work in gathering icing data for use in the prediction of dangerous icing events [7]. Icing data were previously gathered from around 7000 questionnaire responses across countries but were not publicly available [7,11]. Samuelsen et al. report a few cases from the past where the icing data has been made available. Precise icing rate estimation ideally requires detailed information on metocean parameters and ship characteristics during the icing events, which is lacking in existing datasets [7]. This makes the existing full-scale data not detailed enough for validation of computational results [10].

All existing models use observations or empirical equations to estimate several inputs required for predicting sea spray icing. Existing models rely on just three studies, each on one particular ship [7]. Estimation of spray flux is the weakest link in existing models; for example, the spray flux used in ICEMOD is 10–1000 times less than the spray flux used for RIGICE04, and it is difficult to say which one is better [1]. Spray duration and spray period are also derived from limited observations and different researchers have used different formulations for estimating these [7]. Samuelsen et al. compare predictions using two different formulations of spray flux. Other researchers, like Kulyakhtin, have arbitrarily chosen one of the formulations and suggested that any formulation could be used. The performance of individual empirical formulations is beyond the scope of the current study. Details of some of the existing models can be found in Appendix  A.

Due to different approaches, predictions from existing models show considerable variation [10]. Moreover, the limited comparisons with experimental or full-scale data make it difficult to say which model is the most accurate [12]. ISO 35106 notes that none of the existing models predict sea spray icing on a wide range of vessels or structures and that the Overland model has been traditionally used to estimate icing on small- to medium-sized fishing vessels [12].

Detailed full-scale measurements of sea spray icing are expensive [1,10]. Only a handful of these exist and even fewer as published research [1]. Although some researchers like Samuelsen et al. [7] attempted to validate their models against these datasets; the current study found the datasets not feasible for the purpose of validation. This was also observed by Bodaghkhani et al., who pointed out that the spray parameter measurements are not detailed enough for the validation of computational results [13]. One of the reasons is that the icing is heavily dependent on spray flux, i.e., the amount of water impinging on the freezing surface, which varies greatly with the measurement location and could be vastly different at different locations on the same vessel or structure for the same spray event. Taking icing measurements at a different location than the spray measurements would not suffice for validation. A detailed explanation as to why existing datasets could not be used for precise validation is provided in the further sections.

The next best option is to have validation data from controlled experiments—which have also proved to be time-intensive [14]. Owing to the lack of detailed data for sea spray icing, Deshpande et al. carried out controlled experiments with sea spray icing including 20 tests structured into eight experiments. Each experiment focused on one of eight crucial variables affecting sea spray icing. The experiments simulated sea spray icing conditions in a cold climate laboratory. Seawater from a cooling chamber at 2–6 °C was sprayed onto a plate in a freezing chamber at −5 °C to −15 °C with a nozzle. Wind was simulated with a fan in the freezing chamber. Effects on icing rates by varying air and seawater temperatures, wind speed, material, spray duration and period, salinity, and flux were examined. The weight of the ice accreted on the plate was monitored per second for each hour-long test. For complete details, interested readers are encouraged to refer to Deshpande et al. Data for each experiment was presented. Due to the vast nature of the data gathered, the current study explores the use of new techniques including machine learning (ML) to analyze and make sense of the data from the experiments. Adding to the aforementioned challenges in full-scale icing measurement, gathering comparable and varied full-scale data is difficult, especially as each variable necessitates measurements under diverse weather conditions. Scaling effects from lab experiments must ideally be confirmed via full-scale tests. Future research could explore suitable full-scale measurement techniques.

Applications of machine learning algorithms have grown exponentially in the past few years. The technology makes it possible to make sense of large amounts of data that at first look could fail to show any logical relationships. Machine learning has been recently used in icing severity predictions in the aerospace industry [15,16]. Icing on aircraft blades is purely atmospheric icing where the physical icing process, variables involved, and icing severity are completely different from that of sea spray icing. To the best of the author's knowledge, this is the first study wherein a data-driven machine learning model is implemented for the purpose of sea spray icing predictions.

The current study proposes a sea spray icing prediction method based on the experimental data by Deshpande et al. using machine learning and feature engineering techniques and is an extension of the previous study. Post-publication, ten additional repetitive tests were conducted. First, the data are analyzed, followed by implementing some machine learning models without feature engineering, comparing predictions on unseen data. The advantages and disadvantages of these models are discussed in short. Finally, a model, hereafter called the “SPICE” model (prediction model for Sea sPray ICE), is described and implemented with more complexities added in terms of feature engineering in the chosen machine learning model. SPICE is proposed to serve as a new model for the prediction of sea spray icing based on experimentally validated data. Predictions from SPICE are compared to existing models for nine hypothetical cases to see the effects of certain variables on the predictions. Attempts to validate the model based on two full-scale datasets proved insufficient as the datasets were found to be inadequately detailed. Advantages, limitations, and scope for further improvement of SPICE are discussed. The theory behind individual machine learning models is beyond the scope of the current study, but some elementary details are provided for readers not familiar with the subject. Problems with the use of current spray flux formulations are discussed, and suggestions for future research are presented. The SPICE model is, after further improvements, intended to be deployed as a user-friendly application for general use for icing predictions and as a python package for use in further research.

The SPICE model comprises four python classes: data processor, transformer, predictor, and feature calculator, along with their methods. An instance of the data processor manages initial raw data processing. An instance of the transformer computes and integrates new features during feature engineering before the prediction model. The predictor class is a container for the trained ML object and is used for predictions and post-processing. The feature calculator is an addition for cases where the original variables required for predictions are unavailable. The usage of all classes is elaborated in appropriate sections.

Assembling together the methods developed for data processing, feature engineering, the machine learning model, and the feature calculator in Secs. 3–6, respectively, the software architecture for the SPICE model was developed for which the flowcharts are presented in Fig. 1 and Appendix  D. The model is a combination of four different applications:

1. Spice Data Processor: This is used for preprocessing raw data and arranging the data in a format that can be used further by the machine learning model. Details about the data processor are covered in Sec. 3.

2. Spice Transformer: This is used for Feature Engineering. The Spice transformer generates additional data (transforms the dataset) for the ML model. It generates several new columns as a combination of existing input and is used in training and predictions. The need for the Spice Transformer has a basis in the nature of the experimental data and is explained in Sec. 3.4. Details of the Spice Transformer are covered in Sec. 4.

3. Spice Predictor: This is the core ML part of the application which is used while training the ML model and for making predictions from new and unseen input data. While Sec. 3.4 discusses the preliminary ML model evaluation, details of the Spice Predictor are covered in Sec. 5.

4. Spice Feature Calculator: Any ML model requires identical input variables for prediction during training. In Spice, certain prediction inputs (e.g., spray period) can't be directly measured with onboard instrumentation, necessitating estimation from related factors like wavelength, ship speed, and heading using existing formulations from literature. While making predictions, if obligatory variables are absent from the input, the SPICE Feature Calculator is called to estimate them from available data. The flowchart for the Spice Feature Calculator is shown in Fig. 12, and its details are covered in Sec. 6.

Figure 1 shows the flowchart of the SPICE algorithm.

This section gives details about the data used in the current study and how it was processed.

The experimental setup consisted of two climate-controlled chambers for simulating sea spray icing conditions. Seawater from the cooling chamber was sprayed onto a plate in the freezing chamber, on which the ice accretion was measured, with a nozzle. The spray duration and period were controlled with a Programmable Logic Control (PLC) valve. Spray flux was varied with the pump pressure, but also varied due to other factors such as spray duration and spray period. Salinity was varied using different batches of seawater and fresh water. The temperatures in the climate chambers were varied and controlled with the refrigeration mechanism. A fan in the freezing chamber simulated variable wind speeds with a regulator. The detailed setup, procedure, and further details of the experiments and the data collected can be found in Deshpande et al.

Data was recorded every second for each hour-long test. The parameters that were controlled, measured, and used in the study; i.e., “features” for the ML models are thermal conductivity of the substrate material, air temperature, seawater temperature, salinity of seawater, spray flux, spray duration, spray period, and wind speed. The weight of accreted ice was measured and the instantaneous icing rate, also the “target” for the ML models, was calculated from 10-min averaged values of the weight of ice accreted.

The refrigeration mechanism introduced air and seawater temperature variations in each test, acting as covariates across all experiments. This hindered drawing conclusions on the effects of some of the tested variables on the icing rate directly from the experimental data. Other variables remained constant per test, enabling the inclusion of data from every second. It was already proven that the effect of the material diminishes after the initial few minutes (15 min kept as reference) [14], allowing the removal of “material” from the list of features if the initial few minutes are neglected. Data from two “bad” experiments [14] along with the initial 15 min of each test were excluded. Additionally, the last 15 min of each test was omitted to prevent errors in the 10-minute moving average icing rate calculation.

Figure 2 shows the correlations between the variables in the data. Notably, higher wind speeds strongly correlate with higher icing rates with a correlation coefficient of 0.9, while salinity has a minimal impact. Contrary to theoretical expectations [4], experimental data show a slight positive correlation between icing rate and salinity, observed previously by Ryerson [19]. Deshpande et al. acknowledge potential variability even under controlled conditions (Experiment 0). Such variations might contribute to the positive salinity-icing rate correlation, suggesting that errors in salinity would not significantly alter predictions in icing rates. Negative correlations appear for air and water temperatures, and spray period, indicating higher icing rates with lower values. Conversely, positive correlations are evident for wind speed, flux, and spray duration. Some features exhibit a high multicollinearity; e.g., flux shows a high correlation with spray duration and period as these variables directly affect the amount of spray hitting the plate. In a real case, the flux would be also a function of the position on the ship. Air and water temperatures show a relatively high correlation. This could also be true in a real case where the near-surface seawater temperatures vary with the air temperature. The relatively high correlation between wind speed and air temperature can be ignored as they are completely independent variables.

From the correlation analysis, it could be argued to be acceptable to remove salinity, and the spray period forms the list of features due to the low correlation of salinity with the target variable, icing rate, and the high correlation of the spray period to flux. Exclusion of these features was investigated, and the features were retained for better accuracy (more about this in Sec. 5.2).

Figure 3 shows the pair plot with 100 random samples from the dataset for the target variable instantaneous icing rate (instIcingRate). This gives clarity on how the icing rate is affected by each variable and the extremities of feature values. Notably, air and water temperatures spread quite well, while other features are heavily skewed. This is because features other than the temperatures remained constant for an individual test. More details can be found in Table 1 of Deshpande et al. It is later seen in Sec. 3.4 why feature engineering techniques needed to be implemented due to the skewed data.

The dataset contained no missing values and was shuffled and split randomly into train and test sets (X_train: all features, y_train: values of instIcingRate for corresponding X_train rows, X_test: all features, y_test: values of instIcingRate for corresponding X_test rows) where 80% of the rows were used for training and 20% for testing with instIcingRate as the target variable and other variables shown in Fig. 3 as features. The training dataset is used for training the ML models, whereas the test dataset is used for evaluating the performance of the models. Variable units do not affect predictions, but consistency in units is vital during training and predictions.

Table 1 shows the range of all the features in the training dataset. This is very important, as any predictions made for values outside this range would not be said to be accurate.

Practical usage of the SpiceDataProcesser Class is given in Appendix  E.

Four ML models were assessed using the most commonly used metrics [20]: MAE (Mean Absolute Error: best fit value = 0), MSE (Mean Squared Error: best fit value = 0), and R² (R square: best fit value = 1) for the train and test datasets are shown in Table 2. The linear regressor model and three decision tree-based models: (DTR) Decision Tree Regressor, LGBM (Light Gradient Boosting Method) Regressor, and XGB (Extreme Gradient Boosting) Regressor [21] were implemented with the help of the Scikit-learn library [22] with default hyperparameters. The scores in Table 2 are obtained by comparing the predicted instIcingRate values with actual values (y_train or y_test) using the metrics module in Scikit-learn.

Tree-based regression models are preferable over linear regression for skewed and multicollinear datasets [23]. Table 2 confirms this where all metrics show better fits with tree-based models compared to Linear Regression. Thus, linear regression, and additionally DTR, for which the training set shows clear overfitting, are ruled out. XGB Regressor performs best with the lowest MAE and MSE for both sets, serving as a reference (Model 1) for further development. Internal algorithm specifics of individual ML models are beyond the scope of the current study.

XGB Regressor's strong scores in Table 2 indicate satisfactory predictions for unseen data. However, guidelines for evaluating models (metric thresholds) are arbitrary and may not be appropriate in practice [24]. A limitation of this model was identified when testing for u_mod values from 0 to 20, which can be seen in Fig. 6.

Despite good scores, Model 1 showed an unnatural stepped trend in icing rates with wind speed. This is unlikely in reality. Upon examining a tree from the model clarified the reason. Tree-based models categorize data based on rules applied to the feature values, leading to stepped outputs. The severity of this tendency depends on training data, and in this case, has limited practical value. Wind speeds varied from 0–6 m/s in steps of 2, and these were constant throughout an individual test. Transforming wind speed to u_mod created wider intervals. The model broadly categorized data based on u_mod values in the training dataset. The model's categorizations along the average of successive values of u_mod were sufficient for the training data. With consistent u_mod values in the test dataset, the model makes predictions based on these broad categories, yielding good metrics. However, the model underperforms on dissimilar validation data. For example, the wind speed values present in both, the training and testing dataset, were 0,2,4,6 m/s. The model might end up giving similar results if the wind speed in the validation data is 2.5 or 3.5 m/s, which could be practically untrue. The case with wind speeds is especially important looking at the mean absolute SHAP [25] values in Fig. 4. SHAP (SHapely Additive exPlanations) values show feature influence on predictions, but do not infer any causality [26]. Windspeed (u_mod) has the highest impact on the prediction, accentuating stepped appearance versus the icing rates. Other variables might show similar trends but to a lesser extent.

To address predictions between experimental intervals for dataset variables involved exploring additional methods, including feature engineering. Feature engineering encompasses transforming, generating, selecting, or mapping new features from existing variables [27]. An example is when windSpeed was transformed to (u_mod) using domain knowledge. Two algorithms were employed to create new features: RuleFit [28] and FFX (Fast Function eXtraction) symbolic regression [29]. The methodology and positive effects on model predictions using feature engineering with these two algorithms for regression tasks were demonstrated by Cote 2022 [30].

RuleFit, using a tree-based base ML model, generates decision rules as polynomial-like equations from decision paths [31] comprising combinations of original features. Rules are given coefficient and importance values and selected rules are used as new features while training. The python package for the RuleFit algorithm by Molnar [32] was used in this study. FFX produces a single equation that best fits the dataset. This equation adds another feature to the dataset. The python package for FFX by McConaghy [33] was used in this study. Hyperparameter optimization is used for tuning the learning process. A range of hyperparameters is defined, yielding optimal values for training the ML model. This was done using Hyperopt library in python [34]. Grid search is an alternative approach.

RuleFit permits either Gradient Boosting Regressor or Random Forest Regressor (RFR) as the base model. Here, RFR from scikit-learn library [22] was chosen. Hyperopt was used to get optimized RFR hyperparameters. RuleFit was then executed on training data, producing over 100 rules with varying importance and coefficients. Rules with absolute coefficients and importance over a threshold of 0.02 were selected, resulting in 9 rules. Rules with single variables marked “linear” were excluded, as they duplicated original variables. Coefficient and importance values are only used for rule selection. Rules were translated into equations for coding, where variables were multiplied. After coding for clarity and interpretability in python, equations (Table 3) were derived. These equations compute new features from X_train variables, temporarily stored for use in FFX.

Duplicate equations were intended to be removed but were erroneously allowed and discovered at a later stage. However, the power of machine learning models could be seen at a later stage in Fig. 5, where the new model does not consider the duplicates as important and does not severely affect the predictions. In any case, duplicates should ideally be removed in future versions.

The dataframe with existing and new features after RuleFit is used along with y_train to fit the FFX model. The FFX model score was 0.9888. The output of the FFX model is a rather long equation that was cleaned for its use in python and a new feature was calculated. The resulting dataframe is the transformed dataframe that will be used for training the new model. The head of X_train_transformed is shown in Appendix  B for better understanding. Practical usage of the SpiceTransformer Class is given in Appendix  E.

Table 2 shows that the XGB Regressor model gave the best results during the preliminary ML model evaluation. The same model is selected for further use. The only difference is that instead of training with X_train, the model is now trained with the X_train_transformed obtained after feature engineering. Additionally, Hyperopt is used for hyperparameter optimization similar to what was done for the RFR during the training of the transformer model. Practical usage of the SpicePredictor Class is given in Appendix  E.

Table 4 compares the metrics of the model using SPICE (hereafter called Model 2) to Model1. The scores for Model2 show a 17.13% increase in the MSE for the test set. Using SPICE shows improvement in the predictions over XGB with default parameters without feature engineering.

The mean absolute SHAP values for Model 2 are shown in Fig. 5. The new calculated features (see Table 3 and Fig. 10) are given more importance than the original features. The highest mean SHAP values are seen for Eq. (5) feature which is a multiple of air temperature, wind speed, and flux. These are in fact the variables that are most crucial for the prediction of icing rates.

Figure 6 compares the predictions of Model 2 (using SPICE) to Model 1 for the sensitivity of the models to windspeed (u_mod) with constant values for the other features. The implemented feature engineering techniques make the model better at predicting values in between the experimented values of features, in this case, the wind speed. Since new features from the transformed dataframe are used for training, the dependence of the model on individual original features is reduced.

Similar models with the same feature engineering techniques were made to see the effect of the exclusion of certain features that were deemed to be less important by the ML model. The metrics for the different cases are presented in Table 5, where 1 represents inclusion of the feature and 0 represents exclusion.

It was seen that the first case presented in Table 5 with all features included, also the previously named Model 2, showed the best metrics. Thus, all the original features, even though of less importance to the model, were retained for better accuracy. The practical significance of this is that variables like salinity, if unknown, are set to a certain default value, and this assumption would not have any severe consequences on the prediction; however, the prediction would be marginally better if the accurate values of these variables are known.

Until this stage, the model was tested only on arbitrary values of the original features. However, in a real case, many of these values are not straightforward to know. Air and seawater temperatures, and salinity, can be easily measured. As mentioned in Sec. 3.1, the feature u_mod is equated with the relative wind speed, which might or might not be directly available depending on the actual case. Spray duration, period, and flux are reliant on complicated ship-ocean dynamics and may not be directly available and current models use some formulations to estimate these. Spice is provided with a feature calculator to estimate the unknown values before the predictor is called. However, a case for minimizing and as far as possible, not using the feature calculator, especially for the estimation of the spray flux is later argued for. Practical usage of the SpiceFeatureCalculator Class is given in Appendix  E.

Formulations for the estimation of some variables are debated in the literature and are discussed in the following sections. This study does not vouch for the use of any of these formulations and thus it is possible to use any of the available formulations by setting the choice as an attribute of the object (see Fig. 1) of the feature calculator. These attributes are used wherever necessary inside individual feature calculator modules. Each attribute has an arbitrarily chosen default value.

The Horjen formulation (default) or the Beaufort scale is available in SPICE and can be selected for wave height estimation. In case of the Beaufort scale, the mean of the range of wave heights for the corresponding u₁₀ is selected.

Air temperature is an obligatory input to the feature calculator. The feature calculator checks the input data for airTemp and stores it in the temporary dataframe. If airTemp is missing, an error is raised with a notification.

If the water temperature is available, it is stored in the intermediate dataframe. If missing, it is set to 2.14 which is the mean of the values the model was trained on, and the user is notified of the assumption. The water temperature was not set as obligatory since it was later concluded that variation in water temperature caused a negligible change in the icing rates (see Fig. 7).

As discussed in previous sections u_mod is set equal to the relative wind speed (V_r) for a real case. If the relative wind speed is available, this is stored as u_mod in the intermediate dataframe. In the absence of knowledge of the relative wind speed, it is estimated through the submodules. Figure 12 can be referred to for the detailed flowchart.

If available, the spray period is saved in the intermediate dataframe. If absent, the feature calculator checks for the number of sprays per hour (spray frequency) and converts it to spray period or vice versa. If both are missing, submodules assess available data to estimate the time between successive ship-wave collisions (t_col).

The Zakrzewski, Lozowski, and Samuelsen formulations are available in SPICE and could be used by selecting the appropriate attribute. Any divisions by 0 are set to 0.

Various studies employ Liquid Water Content (LWC) in the wave impact spray for flux estimation [3,5,46,47]. Forest et al. noted vast spray flux variations due to different LWC formulations, ranging from 6.94 × 10⁻⁷ to 44 kg/m²/h at 10 m above the mean sea level and attributed this to the difference in climatic factors and spray-generation mechanisms used while developing the LWC models [47]. Kulyakhtin [1] highlighted flux uncertainty as a key icing prediction challenge, pointing out that the flux formulation in ICEMOD gives 10–1000 times less flux than the one in RIGICE04. It is therefore recommended to calculate or measure the flux on a case-to-case basis. Figures 4, 5, and Table 3 show the critical role of flux in icing predictions. Models reliant on formulations for the estimation of spray flux limit these models only to the type of ship or structure the flux formulation is based on.

Nine cases are presented for the validation of the Spice model by comparing the results with Spice to that of five existing models: Overland model (ov), Roebber and Mitten model (rm), (modified) Stallabrass model (st_mod), (modified) MARICE model (ku_mod), and the RIGICE04 model (ri). Details of the existing models and their equations can be referred to in Appendix  A. Each case is set up such that only one variable is varied in the individual case while other variables are constant for that case. For cases 1–7, all features were set within the trained range of the model as given in Table 1 (except for the variables that easily default to a certain value, like water temperature and salinity). For cases 8 and 9, the spray period, spray duration, and spray flux are not provided and thus calculated by the feature calculator.

The right segment of Table 6 shows whether icing rates are a function of the tested variable in all the presented models. Diverse model approaches make direct comparison difficult. It was discussed in previous sections how the skewed data and tree-based models lead to a stepped output. This behavior was vastly reduced with feature engineering, though complete removal might require a vast number of additional tests. Hence, Fig. 7 also presents the trend of the predictions with Spice. The purpose of the trend is to give a rough idea of the behavior of the model and a linear fit was considered. The aim eventually however is to ideally have more data so as to completely avoid the stepped output.

Figure 7 reveals considerable disparities in icing rate predictions by different models. Cases 1–3 depict icing rate variations due to three key variables: wind speed, air temperature, and spray flux. All models exhibit increased icing rates with increasing u₁₀ wind speed. Predictions with SPICE lie in the middle of the spectrum of existing models, closest to that of RIGICE04. While most models project zero icing at zero wind speed, R&M and Spice show some minimal icing. For R&M, Eq. (24) yields positive polynomial values at zero predictors, while for Spice, any input flux maintains slight icing. If the flux was unknown, and thereby calculated with the feature calculator, the flux value at zero wind would be zero, giving no icing. It is thus important to know the feature values for an accurate prediction.

Case 2 shows a clear tendency to increase icing rates with decreasing temperatures. However, the dependence of icing on the air temperature varies vastly from model to model. Predictions with Spice lie somewhere in between all models. At higher temperatures, however, Spice shows the highest icing rates.

In Case 3, only Spice predicts increasing icing rates solely with varying spray flux, where spray duration and period are constant. Spice's ability to achieve this is notable, while other models fail to show similar trends. Notably, in the presented models, flux is not considered by the Overland and R&M models, while the rest incorporate various formulations involving wind speed, wave parameters, ship speed, and height from the mean sea level. This limitation could lead to incorrect estimates of icing rates in existing models, particularly in cases where wave parameters remain constant due to steady wind conditions as the theoretically calculated flux at the fore and aft of the vessel would be constant.

Cases 4–7 explore less impactful variables according to the ML model: water temperature, salinity, spray duration, and spray period. Notably, the Overland and R&M models show high sensitivity to water temperature, and despite the Overland model's criticism about this by Makkonen et al. [49], it remains in use [1]. For the other variables, models either neglected them or exhibited minimal changes in icing rates within the tested range for a constant flux. While Spice did show minor variations in predicted icing rates due to these, the absolute change was negligible. The limited range of tested spray duration (max. 1 s compared to 5.57 s observed by Ryerson [19]) underscores the need for more experimental data to confidently predict a wider input range. Similarly, the negligible but opposing trend in icing rates with changing salinity aligns with observations by Ryerson. Importantly, minor measurement errors in these variables have negligible effects on predictions, overshadowed by dominant factors like wind speed, air temperature, and flux.

For cases 8 and 9, only the modified Stallabrass model considers ship speed and heading. RIGICE04 and MARICE are designed for structures, while Overland and R&M, although developed for vessels, disregard ship speed. Spice relies on local wind speed information, specifically u_mod, which is compared to relative wind speed. In these cases, Spice's feature calculator estimated windspeed (u_mod), spray duration, period, and flux. This involved estimating wave parameters from u₁₀, leading to values outside tested ranges for spray-related variables. While an ML model performs well within its training range, extrapolating beyond it lacks certainty, especially when validation data are absent. Section 7.2 covers the inability to use existing datasets for validation. Thus, the certainty of absolute icing rates for these cases cannot be determined, yet their trends for both cases, in terms of the direction of the slopes, match the modified Stallabrass model, albeit with substantially different magnitude for case 9. Higher ship speeds result in elevated icing rates, likely due to increased relative wind speed, while lower rates coincide with specific ship headings, associated with reduced relative wind speed and flux. Notably, Spice generally yields higher icing rate values than the modified Stallabrass model.

The correct estimation of flux is emphasized in the previous sections as well as by researchers in previous studies [1]. The feature calculator in Spice is provided in situations where a rough estimation of icing rates would suffice. However, for accurate predictions, especially where icing predictions are critical for design and safety purposes, the aim should be to avoid using the feature calculator as much as possible, at least for more important features like wind speed and flux. The feature calculator could be used in case of unknown salinity, water temperature, spray duration, and spray period, which showed comparatively negligible effects on the icing rate, with a minor loss of accuracy.

Figure 8 presents an example where the predictions are heavily dependent on the choice of flux formulation. The three different choices of flux show great variation in the icing rates against wind speed, for a constant air temperature of −10 °C and other parameters set to default. Additionally, many of the values for spray duration, spray period, and especially flux, were seen to be out of range of the tested data—thus calling for additional experimental data to increase the range of confident predictions.

Due to complexities and expenses, few comprehensive icing and spray datasets are accessible [1,10] and even fewer are published. Two openly accessible datasets containing extensive icing rates and related metocean and ship data were explored for validation attempts in real scenarios.

Ryerson [19] gives a detailed description of 39 spray impacts and the total ice accreted in an icing event that lasted for 100 h. However, using this data directly for validation poses challenges. Ice measurements are given at four locations, but only one table includes flux observations. Real flux experienced at the four locations would clearly have been different due to observed differences in icing thickness, even at the same height. The measured ice thickness often exceeded model predictions in Fig. 9. Despite lasting only about 2 min cumulatively, the 39 spray impacts span 100 h, and icing measurements cover the entire period rather than after each event. Ryerson reported fluxes ranging from 5.22 × 10⁻⁴ to 18.62 18.62 kg/m² per event [19]. With such high fluxes, there is a possibility that some accreted ice might be removed in a dynamic process involving multiple accretion and ablation episodes during an icing event [19]; which cannot be captured due to the nature of Ryerson's data. Another graph shows very high wind speeds throughout the 100-hour period. Although most models mostly neglect wind-generated spray, it could affect the observations over a 100-hour period where wave impact-generated spray lasted for cumulatively only 2 min. The graph for the main deck level at the starboard bulkhead shows ice thickness decreasing when wind speed reduces and air temperatures fall below freezing, suggesting potential effects due to wind-generated flux.

The closest comparison could thereby only be achieved with other models as shown in Fig. 9. (The consolidated Ryerson data used as inputs for this calculation is presented in Appendix  C.) This graph converts predicted icing rates to thickness and presents accumulated ice after each spray impact, and eventually, the maximum predicted ice thickness with different models for the Ryerson data. However, even among the models, the closest match to Ryerson data is limited to specific observations at the 02-level starboard bulkhead. This comparison does not provide substantial validation insight for any model, including Spice, in real case scenarios.

In another set of full-scale data, Samuelsen et al., for their MINCOG model, present a case with 37 spray events, of which 14 have an icing rate of 0.33 cm/h [7]. However, scrutiny of the metocean data reveals significant variations in other parameters within these 14 cases. Air temperature varies from −20.2 °C to −1.9 °C, water temperature from −1.3 °C to 6.6 °C, u₁₀ from 9.3 m/s to 25.2 m/s, wave height from 1.5 m to 12 m, and ship speed from 1 m/s to 8.6 m/s, among others, for individual cases. It is unclear whether these presented icing rates are observed or calculated by MINCOG. Nevertheless, the exactly same icing rate for cases with such huge variations of other variables is questionable. It could be a possibility that several parameters are presented in the form of mean values over larger periods of time. In either case, a comparison of the predicted icing rates by any model with this dataset, as presented, did not serve the purpose of validation for a real case.

Kulyakhtin names flux estimation as the biggest cause of sea spray icing prediction errors [1]. Figure 8 shows the stark differences in icing rates from three flux estimation formulations. Additionally, Ryerson from their detailed study of flux and icing measurements mentions that the observed flux data were attempted to be written as general equations as a function of other parameters, but most relationships were complex and not significant [19]. This observation clearly discourages the use of any formulations for the estimation of spray flux. Despite this, researchers have been using a handful of available flux formulations. Models using these formulations, often in the form of LWC, would, at most, predict icing rates for similar ships and the exact locations as ones from the flux measurements. Any attempts to generalize these types of ice prediction models over different types and sizes of ships or structures or use them to estimate the distribution of ice would provide inaccurate results. This effectively deems all theoretical models as inaccurate or limited to very specific cases. The problem was visible while attempting to validate models against full-scale data, since the ice thickness and flux were measured at different locations. This calls for a radical change in how future models approach the estimation of spray flux. Computer fluid dynamics (CFD) techniques and real-time measurement of spray flux at the desired locations are suggested for future research.

• Being developed “bottom-up” from data from controlled experiments, predictions with Spice could be said to be highly accurate within the tested range if the input variables are accurately known. This claim about the accuracy can be made owing to the high model scores while predicting icing rates for the completely unseen test data. Dominant variables like air temperature, wind speed, and spray flux significantly affect icing rates, with others having minor impacts. One could argue that some additional variables, like humidity, are not captured and could affect the icing rates in full-scale testing. Kulyakhtin proved that the in-cloud humidity can be neglected as other terms have a much higher effect on the mass of accreted ice [50].

• Spice can be made available as an easy-to-use package in python or even a deployed version which could be integrated with other systems. This makes Spice extremely user-friendly in comparison to existing models.

• Spice has the ability for continuous improvement by retraining the model with data from additional tests, e.g., for increasing the prediction range or including melting rates.

• Spice can be used as a transient model and could be integrated with the ships’ systems to make real-time icing predictions. The problem, however, lies in the measurement of flux. Investigations into whether instruments such as rain gauges could be used at critical points for the real-time estimation of spray flux are suggested for future research.

• Existing models ignore precipitation as a source of flux, and if accurately known, Spice can estimate combined sea spray and rain-induced icing rates, given its broader training on both salt and freshwater experiments.

• Technical limitations with the generation of high windspeeds in the lab necessitated the preliminary feature engineering technique using the u_mod variable. Predicted icing rates using this technique in Spice fall between those of other models, and ideally should be validated in full scale.

• Spice is built “bottom-up” from data from controlled laboratory experiments and is thus experimentally validated in the laboratory. This is already an advantage over existing models; however, full-scale application of the model could be restricted due to pending full-scale validation. This is a topic of ongoing research. Flux estimation using the feature calculator in general remains the biggest uncertainty for accurate predictions. This said, unlike existing models, Spice is developed to be extremely flexible, such that the user is free to use the most appropriate formulation depending on the case or to use input values from any other source or method of choice. In general, the goal should be to make as little use of the feature calculator as possible. Techniques for the reduction of the feature calculator by using Spice in combination with CFD models are part of ongoing research.

• The original MARICE model incorporates CFD to address the problem of unequal flux at different locations at the same height above sea level but uses existing spray-generation formulations and points out that these spray-generation formulations are the weakest link [1]. Spice, however, offers independence in specifying spray flux, making it possible to combine CFD or measured flux data for more accurate predictions tailored to specific vessel configurations or environmental conditions.

The study introduces a novel sea spray icing prediction model utilizing ML, dubbed “Spice”. The model is based on controlled experiments involving 30 sea spray icing cases in a cold climate laboratory. The performance of various open-source ML algorithms was assessed. Given the skewed nature of the experimental data, feature engineering techniques were employed, yielding promising results. Results from Spice were compared with five existing models in nine hypothetical cases, each involving a variation of one parameter. Another case compared results with different flux formulations. The study concludes by addressing challenges associated with using currently available full-scale data for validation purposes. A comprehensive flowchart of Spice and its development is provided for replicability. This is the first study in the field of sea spray icing wherein a data-driven machine learning model is implemented for the purpose of icing predictions.

The Deshpande et al. experiments lacked definitive conclusions on the impact of certain independent variables due to several covariates. Earlier findings established that material does not affect icing rates beyond initial sprays, while icing rates rise with wind speed and decline with temperature [14]. This study concludes that higher spray flux contributes to increased icing rates within the tested range. Furthermore, variations in water temperature, salinity, spray duration, and period exhibit negligible effects on icing rates compared to dominant factors like air temperature, wind speed, and flux.

Some of the findings from this study are presented as follows:

• Air temperature, wind speed, and spray flux stand out as the dominant parameters affecting icing rate predictions. ML methods confirm their combined importance, aligning with logical reasoning; for example, high wind speeds and flux will not cause icing unless the temperature is below the freezing point.

• Existing models employ spray duration and period solely to estimate flux, not as independent variables. Independently, these variables show minimal impact on icing rates for a constant time-averaged flux within the tested range.

• Linear regression gave poor results and the tree-based XGB Regressor predicted the unseen “test data” well. However, the limited number of tests necessitated the use of feature engineering techniques for which RuleFit and FFX gave satisfactory results.

• Existing models predict vastly different icing rates. The varying approaches and different purposes of the models make it difficult for a direct comparison.

• Validation of models with existing full-scale experimental data proved to be challenging. It is difficult to measure flux and icing at a detailed level on a large ship. However, it is essential to have more detailed full-scale tests for the validation of prediction models, in the absence of which data from controlled laboratory experiments is the most reliable for the development and validation of prediction models.

• A “general” model for sea spray icing is one that can predict icing for a wide range of applications (e.g., for both, structures and for all sizes of vessels). Resolving the issue of flux estimation on a case-to-case basis through CFD or real-time methods and feeding this to Spice could ensure its universal applicability making it a potential candidate as a truly “general” model for sea spray icing prediction. This would also ensure prediction of the distribution of ice which is critical for design optimization and research in this direction is already underway.

The author would like to thank Dr. Edward Lozowski, Dr. Tom Forest, and Dr. Anton Kulyakhtin for sharing the RIGICE04 model and kind email exchanges. The author would also like to thank Tom Yang from Einblick for their time and effort in introducing me to possibilities with machine learning in the Einblick software and assisting me with the case in the initial stages. Finally, the author would like to thank Vedangini Sahajpal for her machine learning lessons in python.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

c =

wave phase speed, m/s

d =

water depth, m

e =

saturation vapor pressure at given temp (subscripts a: at temperature T_a; f: at temperature T_f)

g =

acceleration due to gravity (g = 9.81), m/s²

n =

spray generated at every “nth” ship-wave collision, –

z =

height above mean sea level, m

Α =

ship heading with respect to wind (Straight into wind = 0 deg)

M =

spray flux (subscripts refer to the formulation used), kg/m²/h

N =

number of sprays per hour or Spray frequency per hour/hour

S =

salinity of seawater, Ppt

c* =

empirical constant for spray duration, –

z* =

roughness length or roughness coefficient for calculating the wind speeds by the wind power law, m

t_a =

air temperature, °C

t_col =

time between successive ship-wave collisions. s

t_d =

droplet temperature at impact, °C

t_dur =

spray duration, s

t_f =

freezing point of saline water at salinity S ppt, °C

t_w =

seawater temperature, °C

u_mod =

modified wind speed (windspeed_mod) variable used after preliminary feature engineering (Sec. 3.1), m/s

u₁₀ =

wind speed at height of 10 m above sea surface, m/s

H_s =

significant wave height, m

P_s =

spray period, s

T_s =

wave period, s

V_r =

relative wind speed, m/s

V_rw =

ship speed relative to waves, m/s

V_s =

absolute ship speed, m/s

V_sα =

ship speed in wind direction (into wind is +ve), m/s

V_sβ =

ship speed in wave direction (into waves is +ve), m/s

u_z, u_z0 =

wind speed at heights z and z₀ for the calculation of wind speed by wind power law, m/s

IR =

icing rate (subscripts refer to the model used), as per the given equation, but eventually converted to mm/h and kg/m²/h

PR =

predictor for icing rates by overland, m °C/s

β =

ship heading with respect to waves (straight into waves = 0 deg)

Λ =

wavelength, m

ρ_ice =

density of ice (used only for the conversion of icing rates in terms of mass to thickness and set to 900), kg/m³

Existing models employ diverse approaches and input variables for predicting icing rates. They also serve varying purposes such as RIGICE04 being developed specifically for large stationary structures. This diversity hampers straightforward model comparisons. Nonetheless, selected models were approximated, accommodating coarse assumptions wherever necessary for the evaluation of Spice. Conversion of icing rates among the mentioned models was done to mm/hr and kg/m²/h through straightforward calculations using ice density ρ_ice = 900 kg/m³ [43,51].

MARICE [43] is a time-dependent sea spray icing model utilizing CFD. The original implementation of this model was done with the help of a freezing equation for each cell in the CFD simulation. Unlike Overland and Stallabrass models that rely mainly on metocean and structural parameters, direct comparison of Spice to MARICE is challenging due to the scale at which equations are applied. However, coarse assumptions are made for comparison. The original extensive set of equations are not presented here; readers can refer to the study by Kulyakhtin [43]. For comparison, simplifications include a structural diameter of 0.7 m (matching test plate length by Deshpande et al.), a windward stagnation point at 90 deg, the Schmidt number, originally computed with CFD, is set to 200, and a steady-state version of Eq. (18) is used instead of CFD-calculated flux. Removal of accreted ice in the wave-washing zone is neglected. The results presented here are for coarse comparison, not original MARICE results. The resulting model is termed as “modified” Kulyakhtin model (ku_mod).

RIGICE04 was developed for estimating icing on offshore structures [47] and was employed for comparison using its Excel format. To facilitate the comparison, these settings were applied: large vertical column flag set at 1, large column diameter at 0.7 m, large column top height at 15, large column bottom height at 0, and number of large columns at 1. The met data from Sec. 7 for each condition was inputted in the first row, while the second row was set to arbitrary conditions with temperature above freezing point (1 °C) as required by the program. A 1-hour time difference between the rows yielded hourly icing output. Default relative humidity (80%) and atmospheric pressure (1013.25 bar) from the accompanying metocean data file were maintained. Outputs were in kilogram of ice and ice thickness as functions of height above mean sea level, as bottom height was set at 0. Ice thickness was converted to mass/m² using 900 kg/m³ice density. Only the first and last values are calculated in Rigice04, and the rest linearly spaced for simplicity. Icing rates with the RIGICE04 model are denoted as “IR_ri.”

Figure 10 shows the head (first five rows) of the X_train_transformed dataframe for better understanding of how the dataframe looks like after transformation.

Data from multiple tables from Ryerson [19] were used in the calculation of icing rates in Fig. 9. Some data could be used directly, while some had to be converted to the inputs required by Spice. The consolidated data used as inputs to Spice and other models in Fig. 9 are presented in Table 7.

Figure 11 explains the color keys used in the flowcharts for clarity. The flowchart for the Spice Feature Calculator is shown in Fig. 12.

Figure 13 shows the flowchart for the development of the Spice model, including the training, testing, and validation processes. Figure 14 shows the flowchart for the usage of spice as an end-user of the application. The simplicity of the flowchart in Fig. 14 indicates the user-friendliness of the Spice model once it is deployed for use.

The SpiceDataProcesser Class is provided with a method “data_prep” for obtaining the train-test sets from lab data (Fig. 1). Feature and target names along with the train-test split ratio are set as attributes. The data processer class can be used to retrain the model whenever more experimental data is available.

The SpiceTransformer class inherits properties from scikit-learn's BaseEstimator and TransformerMixin classes. An instance of this class is created, specifying coefficient and importance thresholds for rule selection as attributes. The “fit_transform” method trains the transformer object using X_train and y_train, yielding the transformed training dataframe, X_train_transformed. The object is pickled for future use, and retraining can be done by calling the method again. Due to ML's stochastic nature, rules and FFX equations may vary, even with the same dataset. Once the transformer is trained, the “transform” method is called to transform any dataframe with the original features to a dataframe with features including all equations during testing, validation, and general use.

The SpicePredictor class is a container for the ML model for which XGB Regressor was selected after initial assessment. Training is performed by calling the “fit” method with X_train_transformed and y_train. The trained instance can be pickled. Subsequent “fit” calls retrain the model, with potential variations due to ML stochasticity. Predictions are made using the “predict” method (flowchart in Fig. 1). As experiments by Deshpande et al. lacked “no icing” conditions, post-processing predictions are necessary for which the “post_process” method is called. This sets the predicted icing rate to zero if the air temperature exceeds the water freezing point or if spray flux is zero. The condition for zero icing for zero wind was considered, but the possibility of precipitation-induced flux prevents its implementation.

The SPICE feature calculator class utilizes the “calculate_features” method to handle situations where original feature values are absent during prediction, and rough approximations suffice. With the available input data, the feature calculator stores calculations in an intermediate dataframe and eventually returns a dataframe containing estimated original feature values. The “calculate_features” method encompasses multiple functions organized into modules (as depicted in Fig. 12). Table 8 outlines module divisions. Main modules (Modules 1–7), one for each feature, assess feature values in inputs, while 15 sub-modules are called as needed based on available data. Equations for parameter estimation are detailed in subsequent sections. Missing obligatory inputs triggers an error, notifying users of assumptions made to proceed.