Assessment of added resistance estimates based on monitoring data from a fleet of container vessels

A practical estimation methodology of the mean added resistance in irregular waves is shown, and the present paper provides statistical analyses of estimates for ships in actual conditions. The study merges telemetry data of more than 200 in-service container vessels with ocean re-analysis data from ERA5. Theoretical estimates relying on spectral calculations of added resistance are made for both long-and short-crested waves and are based on a combination of a parametric expression for the wave spectrum and a semi-empirical formula for the added resistance transfer function. The theoretical estimates are compared to predictions from an indirect calculation of added resistance relying on shaft power measurements and empirical estimates of the remaining resistance components. Overall, the comparison reveals a bias in bow oblique waves and higher sea states of the spectral estimates as well as the large variance of the empirically derived predictions — particularly in beam-to-following waves. One of the study’s main findings, confirming previous studies but based on a much larger dataset than in earlier similar studies, is that added resistance assessment based on in-service data is complex due to significant associated uncertainties.


Introduction
Seaworthiness is considered the ultimate criterion in ship design.However, standard practice is to optimize the hull shape for calm water conditions with a service allowance of around 15% to account for the added resistance in actual operating conditions, Shigunov (2017).However, the added resistance is not only vital in ship design but even more so in ship operation: The impact of added resistance on ship safety, economy, and overall life cycle is unquestionable.Therefore, vessel performance analysis and voyage optimization need reliable and precise information about this second order force in actual ship operations.In this line, Tsujimoto et al. (2008) conclude that estimates of added resistance due to waves are essential for correctly evaluating a vessel's emissions and thus its carbon footprint.With tightening rules for the efficient design and operation of ships, the relevance of added resistance is amplified.Recently, IMO (2021) stipulated the EEXI (Energy Efficiency eXisting ship Index) and CII (Carbon Intensity Indicator) regulations, which, in fact, may lead to a new regime of involuntary slow steaming due to Shaft Power Limitations (ShaPoLi) in case of non-compliance.In parallel, rising fuel prices and port congestion may be driving factors for voluntary slow steaming.Hence, added resistance could be of increasing relevance in the years to come due to reduced operating speeds.It is stressed that added resistance is impacts.Shigunov (2017) estimates the added power on ships in a seaway by utilizing a Rankine panel method to calculate the transfer functions of wave drift forces/moments in 3 degrees of freedom.In addition, the consideration of propulsive and engine-related characteristics is a crucial aspect of his work.The calculation procedure for a container carrier is verified using noon report data of two sister vessels.Kim et al. (2017) present a study combining CFD (Computational Fluid Dynamics) results of the added resistance transfer function for the S175 container vessel and the estimation of the speed loss in irregular waves.The latter has been accomplished by applying a unidirectional Pierson Moskowitz spectrum.Unlike the already presented simulationbased studies, the increasing availability of in-service vessel monitoring data allows for more in-depth analyses of added resistance in actual conditions.For instance, Lakshmynarayanana and Hudson (2017) derive the calm water shaft power from high-frequency monitoring data of three sister vessels by strictly filtering for environmental conditions.Their analysis of the added power is then carried out using the initially derived calm water power as a baseline and under the same draft conditions.The results reveal considerable uncertainty in the derived added power and a weak dependency on the sea state.Nielsen et al. (2019) utilize auto-logged high-frequency data of a container vessel in worldwide service and determine the added resistance in an indirect approach.Similarly to the former study, the added resistance is considered the surplus to the empirically derived calm water shaft power compared to the actual measured shaft power.In general, reasonable results have been obtained; however, significant variance and possibly erroneous data in following waves were observed.Furthermore, Dalheim and Steen (2020a) evaluate the in-situ added resistance of an offshore supply vessel operating in the North Sea by considering model test data in both calm water and irregular waves as references.The analysis is based on one year of high-frequency data, and the added resistance has been predicted as non-conservative according to simulated data based on strip theory results.Lang and Mao (2020) propose a semi-empirical method for the prediction of the added resistance transfer function in head waves.In parallel, they apply this approach to high-frequency monitoring data of two separate vessels and use the JONSWAP spectrum together with hindcast data from ERA5 for wave modeling.Similar to Dalheim and Steen (2020a), they found an underestimation and subsequently derived a wave height correction factor determined as 3.5 √   .Taskar and Andersen (2021) assess several known procedures for predicting the added resistance transfer function and verify the obtained results in irregular waves with noon reported and high-frequency data of an oil and a gas tanker.The JONSWAP spectrum was used for modeling the wave environment, and bimodal seaways were considered by means of ERA5 sea state parameters.It stood out that significant variance was visible in the residuals -in higher sea states in particular -despite the detailed data pre-processing methodology.
Overall, considerable uncertainty is a common and dominant issue among all the above-mentioned publications.Both large variance and bias are observed making the isolated analysis of in-service added resistance challenging.Prpić-Oršić et al. (2018) report about several possible sources of uncertainty ranging from erroneous/biased sensor data to limited knowledge of the prevailing sea state.Nonetheless, it is emphasized that with the deployment of continuous monitoring infrastructure onboard vessels, the degree of measurement uncertainty (including human error) reduces drastically compared to manually recorded noon reports (Aldous et al., 2015;Christensen et al., 2018).In addition, the quality of met-ocean data, particularly wind and waves, has improved significantly compared to earlier analyses, often relying on visual observations (de Hauteclocque et al., 2020).Despite improving data quality, the isolation of added resistance in in-service data is still subject to substantial uncertainty.

Motivation and objective
Given the presented literature, it stands out that the majority of studies is limited to one vessel class and, in some cases, even to only a few voyages.Hence, a significant contribution of the present paper lies in the paper's consideration of more than two years of continuous highfrequency monitoring data from 32 vessel classes with 228 container vessels.However, no detailed hull shape or propeller information has been made available, which may be seen as a drawback in terms of accuracy, but this is considered the default case in general vessel performance monitoring by the respective industry.The proposed practical estimation methodology is based on a semi-empirical framework, and only information on the main particulars of the vessels and sea state parameters is required.The approach presented by Mittendorf et al. (2022b) is chosen for the determination of the added resistance transfer function, and the modeling of the ambient wave scenario follows from a parametric directional wave spectrum shaped with parameters obtained from re-analysis data for wind and swell using the public domain ERA5 service.
The paper consists of two parts: (1) Initially, the estimation of the theoretical added resistance in short-and long-crested waves is performed in parallel and compared.In general seakeeping calculations, long-crested waves are typically assumed due to reduced complexity and conservative results (Lloyd, 1989).Thus, it will be studied whether the same holds for the higher order force of added resistance.(2) The second part of the paper is about the comparison of spectral results to empirical predictions, which are derived from in-service shaft power measurements.This empirical added resistance is determined by an indirect methodology, as presented in Nielsen et al. (2019), i.e. by resistance decomposition using empirical methods.Due to the large windage area of container vessels, the wind resistance is considered through experimental data from Andersen (2013) according to ISO 15016 (ISO, 2015).In addition, the sea margins of the individual vessel classes will be presented and analyzed to link ship operation to design.In fact, the respective sea margins of theoretical and empirical methods are shown and discussed.In a dedicated section, different uncertainty sources will be pointed out and partly set into the context of ISO 19030 (ISO, 2016).Ultimately, the present study aims to provide a practical estimation methodology of the mean added resistance in irregular waves based on publicly available wave climate data and an established semi-empirical method for predicting the associated quadratic transfer function.

Composition
In Section 2, both the fleet performance data and metocean hindcast data will be introduced and analyzed.The data pre-processing and filtering makes up the majority of this section.Section 3 has its focus on the spectral estimation methodology and the in-direct empirical calculation procedure.The results are presented and discussed in Section 4, while Section 5 provides dedicated remarks about the inherent uncertainty.Finally, Section 6 summarizes the paper and discusses possibilities for extending work.

General overview of data
The extensive database of fleet telemetry data includes 32 vessel classes with 228 vessels in total and has been acquired over two years between January 2019 and March 2021.As a side note, the time of the measurement campaign is characterized by severe disruptions in the global supply chains due to the COVID-19 pandemic.An overview of the fleet's two major characteristics, i.e. length between perpendiculars   and TEU (Twenty-foot Equivalent Unit) capacity, can be seen in Fig. 1.In view of the regression analyses taken from Kristensen (2015), it is observed that the considered vessels resemble the world's fleet satisfactorily for both Panamax and PostPanamax vessels.The main particulars and other characteristics of the individual vessel classes have been compiled from class information and in-house data sheets.On a separate note, the present fleet shows varying geometric characteristics, such as vertical stems, twin-screw arrangements, and mounted cranes for cargo handling.All of this is useful in a potential discussion about the relevance of the present study and the generalization of its results with respect to ships of other fleets.
In view of Fig. 1, it is appreciated that the considered vessel classes are mostly characterized by   larger 200 m, which underlines the importance of added resistance for smaller values of the relative wave length ∕  with  being the (absolute) wave length.Söding and Shigunov (2015) report that moderate sea states, characterized by wave lengths of 50 m to 150 m, are most commonly encountered by merchant ships.Hence, the bulk of experienced relative wave lengths arranges in the interval ∕  ∈ [0.125, 0.75], acknowledging that this regime of wave lengths is of higher complexity for the calculation of added resistance due to the dominant diffraction component.
The GPS position histories of all considered vessels are depicted in Fig. 2 and show a broad spatial coverage.Moreover, the route path projections of all vessels have been processed by the algorithm presented in Ikonomakis et al. (2022) for correcting any erroneous vessel positions.Still, minor inaccuracies are visible in the GPS data, e.g. in the southern part of Africa or Spain.It is stressed that the majority of the fleet operates the Europe-Far East trade via the Suez canal.

Sensor readings and filtering
The individual sensor readings being part of the CAMS 1 (Central Alarming and Monitoring System) have a sample frequency of 1 Hz.The data aggregated into 10-min samples is sent to shore whenever it is possible to establish a satellite connection.The basic sensor readings comprise GPS information, such as Longitude, Latitude, and Speed Over Ground ().Additionally, engine-related data has been acquired, such as shaft power   and shaft revolutions rpm.Interquartile range filtering has been applied to the distribution of   ∕rpm 3 for the removal of potential outliers in the shaft power measurements.Moreover, sea and air temperatures have been measured on a number 1 CAMS is the onboard data acquisition system of Maersk Line.
of vessels, but for consistency, these two sensor readings were disregarded in this work entirely.Speed Through Water (  ) has been obtained by the Doppler velocity log, and it is noted that different types and manufacturers of this sensor have been used for different vessel classes.Ikonomakis et al. (2021) review several known issues of   measurements and establish a filtering threshold in use of  as a reference, basically given as |  − | < 2 kts, in order to account for possible sea currents.The proposed threshold is also applied herein, and samples with  or   below 5kts are filtered out, as this speed regime is considered maneuvering.It is noted that sequences with frozen sensor recordings or with an apparent offset have been filtered out manually in the case of   and   .Moreover, the vessels were equipped with anemometers in the bow to measure relative wind speed  , and direction   .The measured relative wind speed has been adjusted to a reference height of 10 m, and the recommended equations of ISO 19030 have been used for mapping the relative wind speed  , and direction   to the absolute wind speed   and direction .The loading condition, i.e. the draft at bow   and stern   position, were also measured; however, high-frequency draft measurements are known for their reduced data quality under forward speed, as shown by Gunkel et al. (2018).Therefore, draft information from the loading computer is utilized in this work exclusively for robustness.It is noted that the draft data was linearly interpolated to match the 10-min interval of the remaining dataset.Moreover, large trim angles were removed, and a threshold of  =   −   is defined so that only data for which | | < 2.5 m is kept.In Fig. 3, the distribution of the relative mean draft is depicted for all time instances and considering the entire fleet (228 vessels).
As can be inferred from Fig. 3, the actual experienced mean drafts are highly skewed towards the scantling draft, as the relative draft is normalized by the design draft of the vessels.The average Froude number   of the entire fleet was 0.16 in the observation period; however, the average design Froude number is in the range of 0.22.It becomes clear that some of the vessels operate far off their actual design and contract conditions.
The richness of data allows for rigorous filtering thresholds, which may subsequently enhance the quality of results.Crucially, the data is filtered for consideration of steady conditions exclusively.Dalheim and Steen (2020b) present a detailed methodology for detecting unsteady sequences in high-frequency monitoring data based on the windowwise slope of the considered samples.However, in the present case, 10-min samples are available; therefore, the relative variances of   , , rpm, heading, and rudder angle are used to determine steady sailing conditions.The individual thresholds are empirically derived, and in Fig. 4, a     It is noted that MCR refers to the maximum continuous rating of the main engine.

Sea state information
In the final preprocessing step, the filtered data has been merged with wave data, based on ERA5, according to the corrected GPS data and UTC timestamps.It is stressed that samples with unavailable wave data were disregarded; thus, samples with erroneous GPS data (cf.Fig. 2) are dropped.The ERA5 database is part of the Copernicus EU program (Copernicus, 2020) and provides met-ocean data in hourly intervals at a 0.5 deg spacing in both Longitude and Latitude.In principle, the directional wave spectrum is available at every grid point.However, because of the ample required storage space and increased processing time due to the scale of the present dataset, it has been decided to work exclusively with integral wave parameters.Therefore, the directional spectrum (, ) will be reconstructed using a parametric expression for the wave spectrum (at run time).Herein, it is of interest to model bimodal sea states and, thus, parameters of the wind-wave and swell partitions, as well as of the total mean system, were retrieved.Altogether, 10 parameters were downloaded per valid sample: significant wave height  , (total, wind and swell partition), peak period  , (total, wind and swell partition), mean zero up-crossing period   and wave direction   (total, wind and swell partition).It is noted that the wave direction   is converted to the relative wave direction, as seen from the ship, using the following equation:   =   − , where  denotes the ship's heading, which a gyrocompass has logged.It is noted that  0 = 180 deg.defines head waves and 0 deg.represents following waves.For the spatiotemporal interpolation, a bilinear scheme has been used in the case of all 10 parameters.The fundamentals of interpolating sea state parameters along ships' routes are presented in Nielsen (2021).In Fig. 5, the joint distribution of the total   and   extracted from the present dataset is compared to global wave climate statistics.It is noted that the plot in Fig. 5 2022), in which they list the use of weather routing and good seamanship as reasons for the lower experienced sea states.It is also mentioned that both distributions are characterized by different sources and degrees of epistemic (or systematic) uncertainty.Based on the isolines in Fig. 5(a), it is stated that roughly 80% of the entire data has been acquired under conditions with   < 4 m.With that said, it is emphasized that ERA5 data is, indeed, assumed to be the ground truth in this work, acknowledging the fact that re-analysis data is generally an estimate of the encountered sea state.There are plenty of studies in the literature investigating the agreement of ERA5 data to wave buoy data or to other hindcast vendors, e.g. de Hauteclocque et al. ( 2020), or Nielsen et al. (2022).Therefore the validation of ERA5 data and the discussion of accuracy is beyond the scope of this work.Nonetheless, known related uncertainty sources will be presented in the later discussion herein.

Applied methodologies
The longitudinal wave drift force -the added resistance   -is a second order quantity, i.e.   ∝ 2 , where  denotes the wave amplitude.In general, added resistance is considered a function of ship speed  , ship draft  , mean relative wave heading  0 , wave frequency  and wave amplitude  (Strøm-Tejsen et al., 1973).As Tsujimoto et al. (2008) pointed out, the added resistance of slender-type vessels, e.g.container vessels, is of high complexity and variance due to pronounced flare angles, possibly leading to non-linear wave breaking effects, the influence of a protruding bulb, and also due to higher forward speeds.Several uncertainty sources in experiments and numerical modeling are presented in Shigunov et al. (2018) and Park et al. (2015).The wave correction in ISO 15016 suggests two calculation methods for the assessment of the mean added resistance in irregular waves: (1) A spectral calculation method is described, which relies on an estimate of the added resistance transfer function and a representation of the wave environment in terms of a wave energy density spectrum.(2) The application of the empirical formula STAwave-1 (ITTC, 2014), which expresses a dependency on   and the ship's main particulars, is also mentioned in the respective standard.Overall, the formula has a major resemblance to the Kreitner (1939) formula and reflects a similarly narrow definition range in terms of  0 ,  ,   , and most importantly, ship type.Herein, the more general spectral calculation procedure (1) is chosen, and its intricacies are elaborated in the following subsection.

Spectral calculation of added resistance
In most publications about added resistance on ships in service conditions, the prevalent sea state is represented by a unimodal wave spectrum, assuming long-crested waves, when, in fact, real sea states and their representing wave spectra contain several partitions resulting from wind and swell.Moreover, the waves propagate from directions statistically distributed around a mean vector, i.e. waves are shortcrested.That being said, average parameters as such -the relative mean wave direction in particular -have a reduced physical meaning in multimodal seaways.In fact, only 22% of all samples within the used dataset, cf.Section 2, show a sea state that can be represented by a unimodal wave spectrum 2 highlighting the importance of considering multimodal conditions.Therefore, calculating the mean added resistance in short-crested waves, assuming a bimodal wave spectrum, is carried out using Eq. ( 1).The expression is based on the assumption of stationarity and superposition, i.e. linear wave theory in infinite water depth.
Given Eq. ( 1), it is said that the mean added resistance in shortcrested waves R , is computed by integrating the product of the quadratic transfer function   and the directional wave spectrum  w.r.t.wave frequency  and the (relative) wave direction .In contrast, the mean added resistance in long-crested waves R , is calculated with  integrated for all  and the transfer function   is only computed for the relative mean wave direction  0 .Lloyd (1989) states that considering short-crested waves, in general, smooths the extreme variations of the responses.In addition, it leads to an increased roll response in following waves and increased pitch in beam waves.Due to the effect resulting from the consideration of short-crested waves on the first order responses, it is thought that similar observations may show in the case of the second order added resistance.
The semi-empirical formula proposed by Liu and Papanikolaou (2020) has been developed mainly for practical purposes, for instance, when the detailed hull geometry is unavailable, not to mention that the formula requires low computational effort.Its distinct advantage over other approaches, such as the STAwave-2 (ITTC, 2014), is the validity for all wave heading angles and thus the applicability in short-crested waves.Moreover, the underlying approach of Liu and Papanikolaou (2020) underwent a rigorous validation study (Wang et al., 2021) and is recommended by the ITTC Ship Operation at Sea (SOS) specialist committee for the wave correction during sea trials, ITTC (2021).The Liu and Papanikolaou (2020) formulation applies a correction term for wave reflection in short waves (   ) to the added resistance caused by the motions of the vessel, i.e. radiation (   ).The necessary predictors for the semi-empirical approach are the main particulars, i.e.   , ,   ,   , the block coefficient   and the length of run   and length of entrance   .Additionally, the operating conditions are required, i.e. the longitudinal radius of gyration   and the Froude number  .Lastly, the wave length  and wave direction  have to be specified as discrete vectors for computing   (, ).In the present work, the adaptation by Mittendorf et al. (2022b) of the original formula (Liu and Papanikolaou, 2020) is employed.The formula's parameter vector was calibrated for both slender and blunt-type vessels based on particle swarm optimization as well as a database of experimental data comprising 25 different ships and around 1100 samples.The computed non-dimensional   transfer function is presented in Fig. 6 as a function of the relative wave length ∕  and relative wave direction .Note that the vessel class with a capacity of 15,550 TEU is taken as a case study for visualization.
In view of Fig. 6, it is observed that the formula by Mittendorf et al. (2022b) has a broad definition range and that the added resistance is more pronounced in head-to-beam waves.As it has been stated in Section 2, the added resistance in relatively short waves is of particular interest in this work, and in the formulation of Liu and Papanikolaou (2020), the reflection component is calculated by a modified version of Faltinsen's asymptotic formula.The theoretical details of this near-field approach are provided in Faltinsen et al. (1980).Söding and Shigunov (2015) state that Faltinsen's formula provides reasonable agreement in very short waves compared to other (semi-)empirical methods.From the plot Fig. 6, a minor discontinuity in the vicinity of  = 165 deg. is observed, which may be caused by the piece-wise calculation of the reflection component.It is noted that the spectral calculations were performed with 36 discrete wave directions  ∈ [−, ] rad and a frequency discretization of 50 discrete frequencies  ∈ [0.01, 2] rad∕s.The upper cut-off frequency is justifiable by the low pass filtering effect of ships, i.e. they are irresponsive to high-frequency waves.Following the findings of Mittendorf et al. (2022b), there is an evident lack of experimental data beyond the threshold of ∕  < 0.15.For this reason, the transfer function is padded in this particular regime using the value of   (∕  = 0.15).
The directional wave spectrum is calculated using the exact same frequency and direction discretization as the transfer function and is defined as (, ) = () (, ), i.e. a unidirectional parametric wave spectrum  is multiplied by a spreading function .Herein, a JONSWAP (Joint North Sea Wave Observation Project) type spectrum is used for calculating  to be able to capture the ambient wave spectrum also in regions with limited fetch.Thereby it is understood that two separate spectra for wind and swell sea are superimposed for modeling bimodal sea states.In short, the present approach is inspired by the ten-parameter spectrum of Hogben and Cobb (1986), and the used expression is defined in Eq. ( 2).

𝐸(𝜔, 𝛽)
(2) As can be inferred from Eq. ( 2), a cosine-power type spreading function is applied to the underlying JONSWAP type spectrum   .The spreading function acts as a weighting function for preserving the energy of the spectrum and distributing it around a mean direction according to the spreading parameter   .Moreover, it is noted that  , = 2∕ , and that  denotes the Gamma function.The directional spreading parameter   is defined in Eq. ( 3) and is a function of wave frequency  and   , which is 10 in wind and 25 in swell conditions, as recommended by Goda (2000).
for  >  , . (3) The unidirectional JONSWAP type spectrum   is an adaptation of the two parameter spectrum by Pierson and Moskowitz (1964), which is defined for fully developed seas.However, Hasselmann et al. (1973) found that sea states are never really fully developed, but continue to grow due to non-linear wave-wave interaction.For this reason, the peak enhancement factor  has been introduced to the Pierson Moskowitz spectrum for better agreement to measurement data.The following expressions were adopted from the class guidelines as provided in DNV-GL ( 2018) and the JONSWAP type spectrum   is defined in Eq. ( 4).
with   = 1 − 0.287 ln . (4) The normalizing factor   is applied for conservation of energy.The standard JONSWAP spectrum assumes  = 3.3; and turns into the PM spectrum when  = 1.In addition, the spectral width parameter  is part of the exponent, and the underlying Pierson Moskowitz spectrum    is defined in Eq. ( 5).
As can be seen, the Pierson Moskowitz including the JONSWAP spectrum adopt a  −5 decaying behavior in their tail part.The peak enhancement factor  is in most publications a discrete value, e.g. = 3.3.However, in the present work, an expression dependent on the   and   is used and defined in Eq. ( 6), which has been adopted from DNV-GL (2018).It is noted that separate  values are determined for wind and swell partitions.
It stands out that Eq. ( 6) is a piecewise function and follows the definition range of the JONSWAP spectrum, i.e. 3.6 <   ∕ √   < 5.The spectral width parameter is defined in Eq. ( 7) separately for wind and swell partition.For illustrative purposes, the reconstructed directional spectrum, as obtained from Eq. ( 2), is compared to the directional spectrum from ERA5 for a specific space-time point.In the given case, a test site a few kilometers off the coast of the Faroe islands is used, and the result is shown in Fig. 7.The data applies to the time 3rd April 2016 00:00:00 at 60 deg of Latitude and 350 deg of Longitude, both to the exact degree.Hence, the presented spectra are unaffected by possible uncertainty due to interpolation.In general, good agreement in terms of peak frequency and directional spreading is visible in Fig. 7, even though the directional spectrum from ERA5 shows minor asymmetric directional spreading.Modeling this is obviously not possible using the expression presented in Eq. ( 2).In addition, the spectral spreading is slightly larger in Fig. 7(b).However, it is concluded that the presented parametric expressions are capable of approximating the actual ERA5 2D wave spectrum sufficiently.

Indirect calculation of added resistance
The calculation framework adopted herein is taken partly from Nielsen et al. (2019) and builds upon semi-empirical resistance decomposition, i.e. superposition.The added resistance is calculated under steady conditions and derived by subtracting the calm water resistance from the actual measured in-service resistance.Thus, the determination of the empirical mean added resistance R , , i.e. based on measurement data, is summarized in Eq. ( 8).
As can be seen, the actually experienced resistance force is obtained by multiplying the propulsive force     with the propulsive and transmission efficiencies.Following Kristensen and Lützen (2013), we assume   = 0.98 for the mechanical losses in, e.g.bearings.On the other hand, the determination of   =  0     is generally more complex and depends on the sea state.The available data of propulsive coefficients in waves (as a function of   ,   , and  0 ) is very scarce -for slender container vessels in particular.The experimental studies of Saettone et al. (2021) for a container ship and Yu et al. (2021) for an oil tanker are mentioned for the sake of completeness.The fundamental work of Nakamura and Naito (1977) shows model test data of propulsive coefficients of the S175 container vessel in regular head waves and indicates a quadratic decay of   with increasing wave amplitude.Additional contributors to the decrease in propulsive efficiency are cavitation and wake fluctuations in waves.In irregular waves, Taskar et al. (2019) show that the propulsive efficiency may decrease under wave impact using simulation data.Moreover, Shigunov (2017) assumes   and   to be constant, whereas  0 has a dependency on   and .It is, however, it may be sufficient to assume  0 to be a constant as well for the higher forward speed and relatively low sea states.Additional driving factors for reduced propeller efficiency in waves are possible ventilation and an asymmetric wake profile in quartering waves, as shown by Prpić-Oršić and Faltinsen (2012) and Mikkelsen et al. (2022), respectively.Lastly, the open water propeller curves were not available, and for this reason, the ISO 19030 default recommendation of   = 0.7 is adopted herein.
The studied fleet, cf.Section 2, shows a wide variety regarding hull form and outfitting characteristics.For this reason, the Hollenbach (1998) method has been used for the prediction of the calm water resistance   , , which is mainly due to the applicability for twinscrew vessels.Within the Hollenbach approach,   , is calculated by the sum of the frictional resistance   , which is based on the ITTC'57 correlation line, and the residual resistance   , which is determined by a set of empirical formulae.Their parameters were derived from data of resistance and propulsion model tests of 433 displacement type vessels.The experiments were conducted in the Vienna model basin, and the formula was validated using experimental data from the Hamburg Ship Model Basin (HSVA).In Hollenbach (1997), it is stated that the method shows significantly less variance in contrast to other known methods, which stands out, particularly for twin-screw designs.Additional key advantages of the Hollenbach method include an uncertainty-aware estimate and the relatively small number of required inputs.Moreover, it is valid for different draft conditions and considers the entirety of possible appendages, such as shaft bossings or bow thrusters.The code implementation has been made according to Birk (2019).Shigunov (2017) shows that the wind resistance may surpass the added-wave resistance in lower sea states, which is due to the slender hull shape and the large windage area of container vessels.Shielding effects and large separation bubbles caused by accommodation and container stacks amplify the wind resistance, Molland et al. (2011).Herein, the wind resistance is calculated according to ISO 15016 in Eq. ( 9).
It is stressed that   denotes the frontal area, the density of air is   = 1.2 kg/m 3 and that   is taken from Kristensen and Lützen (2013).The wind resistance coefficients   (  ) are extracted from model test results presented in Andersen (2013) for the fully loaded vessel since this is the most conservative approach and also justified by Fig. 3. Lastly, it is stated that the presented approach disregards added resistance due to steering, nor is any effect of fouling of hull/propeller considered.Nonetheless, for the purpose of comparing the spectral approach, represented by Eq. ( 1) to the indirect, empirical prediction in terms of Eq. ( 8), the above-described methodology appears as reasonable, considering the general degree of uncertainty of in-service analyses of ships.

Results
The following section is divided into three parts: (1) In Section 4.1, results of comparisons of the theoretical long-and short-crested wave computations for added resistance are presented.(2) In Section 4.2, comparisons between the short-crested theoretical estimates and corresponding predictions based on shaft power measurements and empirical resistance decomposition are shown.(3) Lastly, in Section 4.3, the sea margins based on the theoretical and empirical added resistance estimates are conveyed.Due to the different sizes of the individual vessel classes, the results are presented in a non-dimensional form throughout this section.As a definition, from this point, 'mean added resistance' is referred to as just 'added resistance' if not otherwise stated.

Comparison of added resistance estimates based on long-and shortcrested waves
In Fig. 8(a), the deviations between long-and short-crested wave calculations of added resistance are shown in dimensionless form by normalizing the values by the empirical calm water resistance   , (Hollenbach, 1998).Specifically, the relative residuals are presented as a function of significant wave height   , and the samples are colored according to the mean wave heading angle  0 .It is stated that only samples from the 15,550 TEU vessel class are shown, and a corresponding plot for the whole fleet is located in Appendix A.1.In the plot in Fig. 8(a), a heteroskedastic behavior of the residuals can be seen, i.e. their variance increases with higher   values, which underlines the non-linear character of added resistance.Additionally, it becomes clear that the magnitude of R , is higher in head and following waves relative to R , , whereas it is vice versa in the regime of bow oblique and beam waves.This particular structure of the residuals is dependent on  0 and is a result of the directional smoothing in shortcrested wave calculations, i.e. 'averaging' over multiple wave headings.Moreover, in a few cases, the residuals exceed 10% of   , , and large scatter manifests in approximately 20% of the residuals, which has been estimated from the isolines.In fact, the four isolines segregate the kernel density estimate into five layers, which applies to all KDE plots in this paper.
It is found that the variance increases not only with   but it is also dependent on the number of modes (peaks) in the assumed wave spectrum representing the ambient sea state.This particular observation is conveyed in Fig. 8(b) in a stacked histogram, which is applicable to the entire dataset.As described earlier, in this study, wave conditions are referred to as unimodal when either   of wind or swell waves are smaller than 0.1 m.Fig. 8(b) shows the relative residuals based on a unimodal wave spectrum (green bars) and a bimodal wave spectrum (orange bars), respectively.In view of Fig. 8(b), it is observed that the residuals between short-and long-crested wave calculations are characterized by more considerable variance in bimodal sea states.Two possible reasons are: (1) In contrast to unimodal seas, the mean wave direction, as a parameter,  0 is physically less meaningful in a bimodal spectrum, as such a spectrum is composed of two partitions, most likely with different directions of the respective wave systems.However, in the case of R , , the added resistance transfer function is calculated for  0 , i.e. the wave energy is concentrated on an 'artificial' mean wave heading, which explains the significant variance in bimodal conditions between long-and short-crested calculations (considering actual wave directions).(2) Swell partitions usually show higher peak periods and, thus, their energy distribution is -at least for the considered ship sizescloser to the resonance region, which in turn increases added resistance and underlines that swell waves cannot be neglected when calculating added resistance in actual conditions.An advantage of considering short-crested waves is the reduced relevance of the discontinuity in bow oblique waves within the utilized semi-empirical method, which stood out in Fig. 6.Similar observations in short-crested waves have been made in Liu et al. (2022) for the original method (Liu and Papanikolaou, 2020).
In Fig. 9(a), the non-dimensional added resistance in short-crested waves is depicted as a function of the mean relative wave direction, and samples are colored according to   .Whereas, in Fig. 9(b), the nondimensional added resistance is presented for   and colored according to the respective  0 value.It is noted that the 15,550 TEU vessel class is used, and the sea state parameters are taken from the ERA5 hindcast data, whereas   is part of the monitoring data.
In view of Fig. 9(a), it is found that R , exhibits clear dependencies on  0 and   .Qualitatively, a similar behavior in magnitude is visible for  0 when compared to regular waves in Fig. 6, as expected.On the other hand, in Fig. 9(b), it becomes clear that the ships encounter relatively short to medium wave lengths exclusively, since   = 10 s corresponds to ∕  = 0.42 (considering infinite water depth).Even though computations in long-crested waves require only 34% of CPU time in direct comparison to the short-crested case, the latter calculation method is chosen for the remainder of this work.These particular computations do not only include the directional smoothing effect but they also distinguish wind and swell waves with their actual direction, which, in principle, yields a more accurate physical representation of prevailing wave conditions.For the sake of completeness, the residual distributions resulting from short-and long-crested wave calculations are presented in box plots for all vessel classes in Appendix A.1.

Comparison of spectral and empirical estimates of added resistance
In the following, data from the 15,550 TEU class with 8 vessels and 12.3 × 10 4 samples in total is considered.In Fig. 10(a), the relative empirical added resistance, resulting from the methodology presented in Section 3.2, together with the theoretical estimate in short-crested waves, are depicted as a function of the significant wave height   .Moreover, the difference ( ) between   and  is shown in a colormap.It is appreciated that the empirical added resistance is subject to pronounced variance or, rather, uncertainty.It is striking that the corresponding values do take not only positive values but also negative ones, which is seen as unphysical in headto-beam waves and underlines the large degree of uncertainty within the empirical predictions.Conversely, the theoretical added resistance shows roughly a quadratic dependency on wave height, as expected.Additionally, an outlier cluster becomes evident in the proximity of   = 0.8-0.9m, possibly caused by frozen sensors.Furthermore, a dependency of the empirical added resistance on  =   −  is appreciated, which appears reasonable and may be caused by three aspects: (1) The effect of sea currents, (2) reduced sensor accuracy, and (3) insufficient filtering for steady conditions, i.e. some samples influenced by acceleration are still included.In addition, a variance decrease with increasing significant wave height (heteroskedasticity) is visible in Fig. 10(a), which is caused by unbalanced data, i.e. epistemic uncertainty.Multiple sources of uncertainty related to the estimation methodologies and ship telemetry data, in general, will be discussed in Section 5.
In Fig. 10(b), the relative residuals between the short-crested and empirical predictions are assessed corresponding to Fig. 10(a).It is stressed that a similar plot for the entire fleet is attached to this paper in Appendix A.2.To show the actual distribution of the residuals, a kernel density estimate is included.As an immediate observation, major uncertainty stands out in 20% of the data (estimated by isolines).Furthermore, a skewness appears within the residuals dependent on  0 , i.e. the residuals tend to be positive in head-to-bow oblique waves, whereas it is vice versa in beam-to-following waves.In fact, a large variance stands out in the following wave regime.The added resistance in following waves is of high complexity and, for this reason, often neglected both experimentally and numerically.Container ships have a pronounced transom stern, which affects the prediction of added resistance in following waves: Transom sterns may lead to flow separation and breaking waves, which are generally challenging to model, especially by potential flow theory, e.g.Faltinsen's asymptotic formula.Additional challenges in following waves are not only the possible non-linear effect of broaching, i.e. an induced yaw moment due to surf-riding but also the ambiguity in wave loads due to the dispersion of waves, i.e. long waves overtake the ship, whereas short waves are overtaken by the ship.It is stated that conservative estimates of added resistance are obtained in higher sea states via the spectral method, but the large variance in the empirical added resistance and the reduced data availability in the regime   > 3 m impede drawing any firm conclusions.However, Shigunov (2017) shows similar results for two container vessels, i.e. with notable variance and possible overprediction of theoretical estimates in higher sea states.
In Fig. 11, box plots present the basic quantities of the underlying residual distributions of the added resistance calculations separately for each vessel class.It is noted that box plots indicate the interquartile range with its body, the median as well as the maximum and minimum of the distribution, i.e. the 0th percentile and 100th percentile, respectively.As in the previous section, the residuals are calculated relative to   , and for ensuring confidentiality, the individual vessel classes are labeled by an ID number and ID14 corresponds to the 15,550 TEU class.
As can be inferred from Fig. 11, the two present estimation methodologies of added resistance yield values of similar magnitude for the majority of vessel classes, indicated by a residual distribution centered around zero.On the other hand, some vessel classes have a notable offset, i.e. the residuals are heavily biased.Several reasons for the bias are plausible, such as a dependency on ship size, slenderness, additional outfitting equipment (e.g.cranes), or propulsion arrangement (i.e. the number of propellers).However, no common reason has been identified among the outlier vessel classes.To add to the list of reasons, initially, it was thought that the poor agreement in the case of the ID6 class was due to relatively low data availability.However, the ID19 class shows a larger share (2.73%) of the total dataset, and the agreement is still considered relatively poor.Moreover, the variance is not constant for all vessel classes, and significant variance in the estimates appears for some ship types (e.g.ID5), but yet again, no common reason can be connected to that.The conclusion of Tsujimoto et al. (2008) is confirmed herein, i.e. obtaining reliable estimates in beam-to-following waves is challenging -especially for fast and slender vessels.However, in contrast, it is concluded that added resistance in stern quartering waves is, in fact, not negligible due to its magnitude and effect on ship safety.Ultimately, it is reported that the obtained results of the  comparison show major variance but a minor bias, albeit of increasing magnitude in more severe sea states.This is also concluded in studies from the presented literature, e.g.Dalheim and Steen (2020a).The most central and noteworthy point from the previous results and discussions is the large variance in comparing theoretical and empirical added resistance.Generally, ship performance assessments, based on in-service data, are typically associated with large uncertainties, and, as mentioned earlier, in Section 5, multiple uncertainty categories will be presented and discussed.

Calculation of the sea margin
The present study shows not only ramifications to vessel performance analysis or sea trial corrections, but it is also believed that analysis of in-service vessel data may impact the design of new vessels for optimizing several characteristics, such as required engine power, for their actual operational profile.Söding and Shigunov (2015) state that the hull shapes of slender and fast ships (e.g.container ships) are commonly not optimized for added resistance, but the calculation of viable engine margins requires the consideration of added resistance(s).Especially during the initial design phases, a sea margin is needed to estimate the required main engine power correctly.As mentioned in Section 1.1, constant values of approximately 15% are most frequently employed for this purpose.Due to the richness of data, it is attempted to postulate general hypotheses for container vessel design and their requirements for engine power in actual conditions.However, seamanship, maintenance scheduling, and the fidelity of the weather routing software are company-specific, which introduces certain biases into the present dataset.In the present case, a dependency of the empirical and theoretical sea margins on the vessels' slenderness ratio is shown in Fig. 12 using the 10th percentile, 50th percentile, and 90th percentile of the individual distributions of all 22 vessel classes.Herein, the sea margin is calculated by dividing the sum of added resistance in short-crested waves and wind resistance by the empirical calm water resistance, i.e.R +    , .As a side note, the excess resistances due to marine growth and shallow water are neglected in the present analysis.
As shown in Fig. 12, both the theoretical-and empirical-based sea margins show dependencies on the slenderness ratio.In general, linear relationships between sea margins and   ∕∇ 1∕3 are visualized using regression lines for all three quantiles .Overall, it appears that theoretically and empirically derived sea margins show an opposing dependency on the slenderness ratio.In general, it is believed that the added resistance, and thus the sea margin, will be reduced for more slender vessels, which can be seen in Fig. 12(a).Moreover, the variance is relatively large in relation to the slope of the regression lines in Fig. 12(b).Hence, no defensible conclusions are possible in the case of the empirically derived sea margin.Moreover, the distribution of the sea margin appears to be asymmetric and skewed towards lower values resulting from the exponential   distribution.The 90th percentiles exhibit significant scatter and samples in both Figs.12(a) and 12(b), which are exceeding a 60% sea margin.Simply speaking, the lump sum corrections recommended for the use in the first design iterations appear reasonable given the average sea margin in Fig. 12 noting that Kristensen and Lützen (2013) suggest sea margins of around 20%-25% for the trade from Southeast Asia to Europe.It goes without saying that in advanced design iterations, more detailed approaches, as presented in Liu et al. (2022) or a simulation-driven design process, as shown in Harries et al. (2019), are needed for viable ship designs.In the end, the analysis of sea margins is of even greater importance in the case of slower and blunt-type ships, such as bulk carriers and tankers, due to the risk of losing maneuverability in severe sea states and thus for the minimum propulsion power requirements.

Discussion of uncertainty sources
Shigunov (2017) states that a large set of influencing factors of added resistance is associated with significant uncertainty or may even be completely unknown.It has been stressed that it was not possible to pinpoint one overall reason for the relatively large uncertainty in the respective estimates of added resistance and their comparison.Hence, several possible contributing reasons are discussed in the following, and several potential sources of uncertainty are presented, with the spectral calculation and the indirect empirical methodology addressed separately.First, however, common influencing factors of uncertainty are pointed out.
The empirical methodology is based on resistance decomposition, and the spectral method relies on the superposition of several wave components.Therefore, both adhere to linear theory and neglect interaction effects per definition.In reality, a ship's behavior in a seaway is non-linear, especially in higher sea states, and, therefore, superposition, in general, is not strictly valid.As a matter of fact, non-linear phenomena, such as propeller load fluctuations, propeller ventilation, and slamming (resulting in whipping-induced vibrations), increase in likelihood and magnitude in higher sea states.Moreover, the wave-induced transverse drift force and the yaw moment have been disregarded in the present work, but Shigunov (2017) reported a minor effect of an unrestrained heading on the results of added power in waves.The chosen physical parameters may also introduce uncertainty in both methods.For instance, seawater properties for a temperature of 15 • C, i.e.   = 1025.9kg/m 3 and kinematic viscosity  = 1.19298 × 10 −6 m 2 /s, were assumed to be constant, according to ISO 19030, when in fact both are dependent on temperature and salinity.Therefore, it is believed that the viscous resistance component is subject to scatter and possibly overestimated in regions with higher temperatures.Furthermore, in this study, all ships are represented by only their basic ship main particulars; since no details about hull or propeller shape or exact loading conditions were accessible.

Uncertainty sources within the spectral method
Inherently, the spectral calculations of added resistance do not capture the interaction effects and unsteady phenomena due to the underlying assumptions of Eq. (1).Moreover, linear wave theory assumes moderate wave steepness and it has been shown experimentally in regular waves in Park et al. (2019) that wave steepness has a significant effect on the added resistance transfer function in short wavelengths and on the occurrence of wave breaking.Therefore, the effect of wave steepness,   = 2   2

𝑧
, in irregular waves on the residuals has been studied in Fig. 13.
In fact, Fig. 13 reveals no significant correlation between residuals and steepness, but the variance of the relative residuals increases notably with higher   values.A clear error source within this analysis is that a mean steepness aggregating both wind and swell partition is used.Furthermore, the application of simplified expressions for both  the quadratic   transfer function and the parametric wave spectrum come with their own uncertainties.For instance, the present parametric directional spectrum is limited to symmetric directional spreading, while directional wave spectra of natural ocean waves rarely appear symmetric in their energy density directional distribution (Barstow et al., 2005).The retrieved sea state data and its interpolation introduce additional variance, as discussed by Nielsen (2021).Moreover, de Hauteclocque et al. ( 2020) compare several hindcast vendors (including ERA5) to data from 9 wave buoys and show that ERA5 data tends to underestimate the wave height by approximately 4%.In the same light, the comparison of ERA5 to altimeter (satellite) data reveals a reduced underestimation of only 1.6%; however, more significant deviations have been identified in harsh conditions, i.e.   > 8 m.In the present context, such conditions are mostly circumvented by weather routing; therefore, the associated discrepancies are not necessarily influential.In Section 4.1, different behavior in the residual variance has been observed in uni-and bimodal sea states, and, in fact, similar findings as in Fig. 8(b), but with larger variance, have been made in case of the comparison of theoretical and empirical added resistance.The corresponding histograms are placed in Appendix A.2.The employed semi-empirical method of Mittendorf et al. (2022b) is already characterized by a reduced model parameter uncertainty due to the conducted calibration procedure based on particle swarm optimization and experimental data of slender and blunt-type ships.Still, notable variance and underestimation were observed in Mittendorf et al. (2022b) while validating the method for a fast and slender ship.As reported by Tsujimoto et al. (2008), the estimation of added resistance of slender vessels is problematic due to their pronounced flare angles (possibly leading to wave breaking), their higher operating speeds, and their distinct bulbous bow.Furthermore, the exact loading conditions and, thus, the metacentric height and the longitudinal radius of gyration   are unknown.As shown by Holt and Nielsen (2021), the latter has a sizable impact on the added resistance transfer function and has been approximated with   = 0.25  throughout this study.In fact, Dalheim and Steen (2020a) conclude that both limited knowledge of the actual   and nonlinear effects, such as wave breaking, are leading contributors to the deviations in their comparison of theoretical and empirical added resistance.

Uncertainty sources within the empirical method
The indirect empirical predictions are based on superimposing several resistance components, and the underlying calm water resistance method relies on a similarity to the parent hulls of the Hollenbach (1998) method.However, specific geometric characteristics, such as vertical bow shapes, as in the case of the ID1 class, were not considered among these underlying ship models.In the future, it could be interesting to study the sensitivity to several different empirical calm water resistance methods, such as the approach by Guldhammer and Harvald (1974).Furthermore, changes in engine characteristics in waves were neglected throughout this study, even though transient behavior, such as ventilation and propeller racing, may lead to engine overload, as described in Prpić-Oršić and Faltinsen (2012).Clearly, calculating the propulsive efficiency   is essential for the empirical calculation methodology.Even though Shigunov (2017) found that the influence of   on actual engine power in a seaway is one order of magnitude less than the influence of a resistance increase due to waves, which in turn justifies, to some extent, that the current focus is merely attributed to the accurate modeling of the wave-induced added resistance.Still, it has been shown by Nakamura and Naito (1977) that   decreases quadratically with increasing   , and consequently, the ISO 19030 default value of 0.7 loses its validity in higher sea states.Therefore, the found deviations in higher sea states in Fig. 10(b) may be partly attributed to the application of the constant value of   .More accurate modeling of propulsive coefficients and propeller characteristics is an essential area for future work.On a separate note, it is thought that data quality of monitoring data is another significant driving factor of the observed uncertainty -despite the rigorous filtering procedure.When it comes to the measurement uncertainty of ship performance data, Aldous et al. (2015) conclude that both speed and power measurements are most affected by measurement error and sensor drift when compared to other essential sensor readings.From Eq. ( 8), it can be inferred that both are vital for the empirical determination of added resistance.Hence, sensor maintenance and quality are of great importance for obtaining sufficient information on added resistance from in-situ data.In Fig. 10(a), we see that the empirical added resistance depends on  =   − .Thus, it is herein attempted to analyze the variance and bias of the different Doppler velocity logs with respect to  as a reference.Obviously, the analysis assumes both quantities to be time-invariant when, in fact, both quantities of   are subject to sensor drift and calibration.However, an estimate of bias and variance may still be indicative of deviations observed in Fig. 11.In Fig. 14, the normalized variance and mean of  are shown for each vessel class.In parallel, the variance and mean (or bias) of the residuals are visualized by the size of the dots and the colormap, respectively.
From Fig. 14, it is appreciated that there is a minor coupling between the variance and bias of  and the residuals of theoretical and empirical added resistance.For instance, the vessel classes of ID6 and ID19 are characterized by large biases in both cases.It is striking that vessel classes with a large bias of the residuals show a negative bias of  .Conversely, vessels with a positive  bias show a negligible bias of the residuals.Broadly speaking, the sizes of the dots (i.e. the residual variance) also increase with higher variance of  .However, ID11 is considered an outlier.All in all, the sensor quality, possible sea currents, and the applied method for filtering for steady samples may be responsible for these observations.It is believed that similar observations can be made for the torsion meter, i.e. the sensor for the brake power   .However, in this case, there is no direct reference value for a similar investigation like for   .Data quality and sensor maintenance are critical issues in the field of performance monitoring.Therefore, it could be an interesting aspect to follow the approach of Sogihara (2021) and study the effect of individual sensor readings as well as wave parameters on the overall uncertainty within the added resistance estimate using a Monte Carlo simulation methodology.
It has been stated in Shigunov (2017) that the wind resistance of container ships can be of higher magnitude than added resistance due to waves in lower sea states, which is due to the large frontal area of container vessels.In fact, ISO 19030 imposes a filtering threshold based on the absolute wind speed   alone for disregarding any instances with more severe environmental conditions, which may lead to higher added resistance due to wind and waves.For fully developed wind seas with unlimited fetch, the relationship   ∝  2  holds, and there are well-known formulae from, e.g.Bretschneider   = 0.0248 2  (Michel, 1999) and   = 0.115 1.41  (Shigunov, 2017), which are visualized in Fig. 15.In the following, the absolute wind speed   , as derived from relative wind measurements on all of the considered vessels, and corrected to the reference height of 10 m, is depicted for the significant wave height   .The samples are colored according to the zero upcrossing period   , and the wind threshold (i.e.7.9 m/s) of ISO 19030 is included as a broken red line.It is noted that the two sea state parameters are subject to potential inaccuracies due to spatiotemporal interpolation and the inherent uncertainty of the ERA5 database.
Generally, it can be seen in Fig. 15 that the kernel density, i.e. the bulk of the data, clusters around the two empirical formulae in the lower sea states, i.e.   < 4 m.Still, significant variance stands out, and even possible erroneous samples are visible, e.g.samples with   > 20 m∕s and   < 1 m.In fact, wind measurements are subject to major uncertainty, which is often a result of the placement of the anemometer and the resulting wake from, e.g.container stacks or the accommodation.Sogihara (2021) reports a measurement accuracy of ±5% and ±5.0 deg. in case of relative wind speed and direction, respectively.In fact, too high wind resistance values were also reported by Schmode et al. (2018) and may be responsible for the negative bias of the residuals seen in Fig. 10(b).In this light, it may also be argued that the ISO 19030 wind threshold could be deemed inappropriate not only because of the weaker relationship of   and   in higher sea states but also due to the fact that swell waves are independent of the prevailing wind speed but yet may be significant in terms of added resistance (indicated by higher   values).
ISO 19030 is heavily discussed in the practical and academic domain and Schmode et al. (2018) suggest several improvements for ISO 19030 and point out drawbacks of the standard procedure, such as the speed and draft dependency of the ISO 19030 performance indicators.Moreover, Schmode et al. (2018) state that the strict filtering mandated by ISO 19030 may lead, in some cases, to the fact that a great chunk of the data has to be disregarded, which directly impacts the reliability of the obtained performance indicators.In view of Figs. 15 and 10(b), it appears favorable to filter samples according to shaft power measurements and a speed-power baseline, e.g. a sea trial curve.Bertram (2016) suggests disregarding any samples with a relative added power (dividing by the reference power from the baseline) that exceeds 100%.
In doing so, not only higher values of added wave and wind resistances (and possibly non-linear behavior) are neglected, but also potential measurement errors are filtered out.Another filtering technique could be based on the estimate of hindcasted   values and thus on the Beaufort scale.In ISO 15016, it is recommended (for   > 100 m) to perform sea trials only if Bft< 6, i.e   < 4 m and   < 13.8 m∕s.A similar threshold appears as reasonable in the case of ship performance analysis, i.e.ISO 19030.
Wind forces as such and their exploitation received increased focus in recent years, both due to the Suez canal blockage in 2021 and due to the re-emergence of wind(-assisted) ship propulsion.Hence, it is believed that dedicated studies for accurate modeling of the wind forces on container (and other) vessels and their uncertainty quantification could be of significant scientific and practical relevance.As pointed out, wind speed and direction measurements may be subject to considerable error due to the distorted airflow around the ship's hull, accommodation, and container stacks.Moreover, major differences in data quality of measurements obtained from different wind anemometers have been found in, e.g., the study of Oh et al. (2018).Therefore, a similar sensitivity study to that of Ikonomakis et al. (2021) for   appears to be highly relevant for ship-based wind measurements and could be carried out using ERA5 wind data as a reference.Such studies may enhance understanding and corresponding estimation procedures of wind loads for both ship performance and sea trial assessment.

Conclusions and future work
This paper assesses theoretical and empirical estimates of added resistance using in-service data from a fleet comprising more than 200 container vessels.The calculation of the theoretical added resistance has been carried out in the spectral domain using an established semiempirical formula for arbitrary wave heading angles (Mittendorf et al., 2022b).The sea state, on the other hand, has been modeled by a parametric expression for a directional sea state based on a JONSWAP type spectrum (Hasselmann et al., 1973), and the required parameters for both wind and swell partitions have been retrieved from the ERA5 database (Hersbach et al., 2018).For comparison, an empirical added resistance prediction was made using resistance decomposition and ship telemetry data, including shaft power measurements.The ship telemetry data underwent a rigorous filtering procedure beforehand to disregard any instances with, e.g.shallow water or unsteady conditions.Similar to other contributions from the literature, pronounced variance, and uncertainty in the comparison of theoretical and empirical added resistance estimates have been found, and multiple uncertainty sources have been discussed.It has been shown that data quality is vital for reliable predictions of   , and the fidelity of both sensors and filtering techniques are essential.Ultimately, as the main finding of the present study, the conclusion of Bertram (2016) is confirmed, i.e. added resistance is generally difficult to predict in actual conditions -particularly in short and oblique waves.
In conclusion, the assessment made in this study calls for an uncertainty-aware estimation methodology of added resistance.In fact, the method of Mittendorf et al. (2022b) for the added resistance transfer function, and the approach by Hollenbach (1998) for the calm water resistance are both capable of providing transparent uncertainty bounds.However, for a consistent estimation procedure, reliable uncertainty information of sea state data is necessary.Following the work of Bos (2018), an ensemble of hindcast providers could form the basis of uncertainty quantification of wave data.Moreover, information on the prevailing wave environment in terms of a wave spectrum with increased accuracy is an essential aspect for calculating the mean added resistance in irregular waves.The application of the wave buoy analogy (e.g., Nielsen, 2017) -as presented in Mittendorf et al. (2022c) -for data acquisition (and for real-time decision support) shows potential for enhancing the quality of sea state information and thereby of in-situ added resistance estimates.On a separate note, it is

Fig. 1 .
Fig. 1.Relationship between ship length (  ) and capacity in terms of TEU considering the individual vessel classes.It is noted that the regression curves are taken fromKristensen (2015) and that  denotes the beam of the vessel.

Fig. 2 .
Fig. 2. GPS position projection of 228 vessels in the time between January 2019 and March 2021.

Fig. 3 .
Fig. 3. Histogram of the relative draft according to the respective design drafts of the vessel classes.

Fig. 4 .
Fig. 4. Filtering methodology applied to one case ship (15,550 TEU) in October 2020.It is noted that MCR refers to the maximum continuous rating of the main engine.
(a) is an aggregated scatter, kernel density mapping, whereas Fig. 5(b) shows a contour plot.The comparison of Figs.5(a) and 5(b) reveals that the observed sea states are skewed towards lower values, notably for both   and   , in contrast to the long-term design probability distribution (DNV-GL, 2018).This finding is in line with the work of Nielsen and Ikonomakis (2021), and Miratsu et al. (

Fig. 5 .
Fig. 5. Joint distribution of   and   (a) according to the present dataset and (b) to long-term statistics.

Fig. 6 .
Fig. 6.The non-dimensional   transfer function for   = 0.16 for one case vessel class (15,550 TEU) according to the formulation of Mittendorf et al. (2022b).

Fig. 7 .
Fig. 7.It is noted that both spectra correspond to the same time and position (  = 3.7 m).Also, it is stressed that  =  2 is displayed and that the spectral density is normalized.

Fig. 8 .
Fig. 8. Relative residuals of added resistance according to short-and long-crested waves.

Fig. 9 .
Fig. 9. Dependency of the non-dimensional added resistance on forward speed (  ), mean relative wave heading  0 (a) and zero up-crossing period   (b).

Fig. 10 .
Fig. 10.Analysis of relative theoretical and empirical added resistance for the 15,550 TEU vessel class in direct comparison (a) and their residuals (b).

Fig. 11 .
Fig. 11.Box plot of the relative residual distribution for all 22 considered vessel classes.Note that the relative samplesize is shown on the left hand side in the plot and that outliers are not depicted.

Fig. 12 .
Fig. 12.Three quantiles  of the sea margin distributions as a function of the slenderness ratio for all 22 considered vessel classes.It is noted that ∇ is the volume displacement at design draft.

Fig. 13 .
Fig. 13.Sensitivity of the relative residuals to the wave steepness irregular waves.

Fig. 14 .
Fig. 14.Normalized variance and bias of  =   − .It is noted that the dots are colored according to the bias of the residuals and their size corresponds to the residuals' variance, as presented in Fig. 11.

Fig. 15 .
Fig. 15.Absolute wind speed   as a function of significant wave height   including two regression analyses and the ISO 19030 filtering threshold (red).

Fig. 19 .
Fig. 19.Note that the ordinate is in logarithmic scale and that bimodal (orange) and unimodal conditions (green) are depicted separately in stacked histograms.