Sea-State-Conditioned Motion Response of Berthed Ships Using Field Measurements from Multiple Vessels and Berths

1. Introduction

Ship motions at berth have long been recognised as a critical factor influencing port operability, mooring integrity, and cargo-handling efficiency. Excessive motion can disrupt terminal operations, increase loads on mooring and fender systems, and compromise safety. In severe cases, adverse environmental conditions have led to temporary suspensions of port activities and associated logistical disruptions. Consequently, the prediction and assessment of berth motions have been widely investigated using analytical, numerical, experimental, and field-based approaches.

Early investigations of ship motions at berth are predominantly based on frequency-domain formulations grounded in linear potential-flow theory. Within this classical linear framework, vessel responses are characterised through Response Amplitude Operators (RAOs), which relate harmonic wave excitation to motion amplitude in each degree of freedom [1,2,3]. For irregular waves, the motion response spectrum $S_{\eta }\left(\omega \right)$ is expressed by Equation (1).

$$S_{\eta }\left(\omega \right)={\left|\mathrm{RAO}\left(\omega \right)\right|}^{2}S\left(\omega \right)$$

(1)

where $S\left(\omega \right)$ denotes the incident wave spectrum, and it comes from linear systems applied to ship motions in waves [1]. The motion variance is obtained through spectral integration of the response spectrum over frequency [1,4], as in Equation (2).

$${\sigma }^2={\int }_0^{\infty }{\left|\mathrm{RAO}\left(\omega \right)\right|}^{2}S\left(\omega \right),\mathrm{d}\omega$$

(2)

By evaluating these expressions for prescribed combinations of significant wave height ${\mathrm{H}}_s$, peak wave period ${\mathrm{T}}_p$, and wave direction, deterministic response envelopes can be constructed for a specified vessel–mooring–bathymetry configuration. Such RAO-based formulations form the theoretical foundation of ship and offshore hydrodynamics [5,6,7]. They provide physically rigorous prediction capability and detailed insight into resonance mechanisms, low-frequency excitation, and hydrodynamic coupling effects.

However, the application of deterministic RAO-based formulations in Operational berth environments present substantial practical challenges. Construction of response envelopes requires detailed knowledge of vessel geometry, mass properties, hydrodynamic coefficients, mooring stiffness, damping characteristics, and local bathymetry. In long-term berth monitoring datasets involving multiple vessels, varying mooring arrangements, and incomplete instrumentation records, such parameters are typically unavailable or inconsistent across observations.

Consequently, a methodological gap arises between physically rigorous hydrodynamic prediction frameworks and the realities of operational berth monitoring. Field datasets are frequently heterogeneous, comprising multiple vessels, varying mooring configurations, incomplete instrumentation records, and non-overlapping measurement periods collected under evolving harbour wave conditions. These characteristics limit the direct applicability of deterministic reconstruction approaches and motivate the development of transparent statistical frameworks that can extract engineering-relevant response characteristics directly from heterogeneous field data.

A further limitation of many conventional frequency-domain and inverse-response approaches is the assumption of statistical stationarity of the wave field over the analysis interval. Stationarity implies that statistical properties such as wave energy, spectral shape, and dominant direction remain approximately constant in time. While this assumption may be reasonable for short-duration records in open-sea environments, it is often violated in harbour and berth settings, where wave conditions can evolve due to varying wind forcing, harbour-geometry effects, reflection, diffraction, vessel traffic, and long-period oscillations.

For this reason, ocean waves are commonly described using statistical characteristics rather than deterministic realisations. As noted by Holthuijsen [8], wave conditions in coastal engineering are typically represented through spectral descriptors such as significant wave height and characteristic wave period. Conditioning berth motions on such Sea-state parameters are therefore consistent with the prevailing statistical wave-modelling paradigm.

Physical model testing has also played a prominent role in investigating berth dynamics. Scale-model experiments conducted in wave basins have been used to study the effects of wave directionality, berth layout, and mooring configuration on ship motions [5,9,10]. While laboratory studies allow for controlled investigation of specific mechanisms, they are inherently limited in representing the full range of environmental variability encountered in real ports, and scale effects may influence the modelling of damping and mooring behaviour.

Field measurements provide a complementary perspective by capturing ship response under real operating conditions, naturally accounting for the combined effects of waves, wind, currents, berth geometry, and mooring systems. Several authors have reported field observations and analysis of ship motions and harbour wave conditions at berth, and compared measured responses with numerical predictions or empirical models [11,12,13,14]. These studies have demonstrated both the value and the challenges of field data, highlighting issues such as sensor noise, data gaps, and variability in environmental forcing. Importantly, field datasets are often heterogeneous, comprising measurements from multiple vessels, different berth locations, and non-overlapping time periods.

A standard limitation of many existing field-based studies is the implicit assumption that multiple DoF motions are available simultaneously for a given vessel and sea state. This assumption enables multivariate analyses and, in some cases, inverse estimation of wave characteristics from ship motions [15,16]. In practice, however, operational monitoring campaigns frequently yield incomplete datasets: different ships may be instrumented differently, sensors may operate intermittently, and measurement periods may not coincide across berths. Enforcing simultaneity across all DoF in such cases can lead to substantial data loss and biased characterisation of the response.

To address these challenges, several authors have advocated statistical or empirical approaches that condition ship response on measured environmental parameters rather than attempting a full deterministic reconstruction [4,12,17]. In these approaches, motion statistics are analysed as functions of sea-state descriptors such as significant wave height, peak wave period, and wave direction. Such methods are particularly well-suited to berth environments, where wave climates are often site-specific and influenced by local sheltering, reflection, and diffraction effects.

Despite the maturity of deterministic hydrodynamic modelling techniques, a methodological gap persists between physically rigorous RAO-based prediction frameworks and the realities of long-term operational berth monitoring. Existing approaches typically require detailed vessel-specific parameters and simultaneous multi-DoF measurements under stationary conditions. In practice, however, operational datasets are frequently heterogeneous, comprising multiple vessels, varying mooring configurations, incomplete instrumentation records, and non-overlapping measurement periods collected under time-varying harbour wave conditions. These characteristics limit the direct applicability of deterministic reconstruction methods and highlight the need for a transparent, statistically grounded framework that can extract engineering-relevant response characterisation directly from heterogeneous field data.

Building on this perspective, the present study proposes a sea-state-conditioned statistical framework for analysing berth-motion data from multiple ships and locations. Each motion degree of freedom is analysed independently, avoiding unnecessary assumptions regarding simultaneity or coupling between responses. A transparent quality control procedure is applied to ensure the plausibility and consistency between observed sea states and measured motions. Motion response envelopes are derived by binning observations in sea-state space and computing representative and conservative statistics.

In addition, synthetic sea surface elevations corresponding to each observed sea state are generated using a spectral random-phase approach [4,8]. These realisations are statistically consistent with the prescribed sea-state parameters and are used to support time-domain interpretation and visualisation of the response envelopes.

Contributions of This Study 

The main contributions of this study can be summarised as follows:

• A data-driven framework for heterogeneous berth motion datasets. A statistical analysis methodology is developed that explicitly accounts for incomplete overlap between motion DoF, enabling the effective use of field data collected from multiple ships and berth locations without enforcing unrealistic simultaneity assumptions.
• A transparent and robust quality control procedure. A two-stage quality control approach combining physical plausibility checks with robust regression-based outlier detection is proposed to identify and remove inconsistent response–sea-state pairs while preserving the dominant field-observed behaviour.
• Sea-state-conditioned motion response envelopes for berthed ships. Motion response envelopes are derived for individual DoF by conditioning observed motion statistics on significant wave height, peak period, and wave direction, providing representative and conservative measures suitable for berth operability assessment.
• Generation of statistically consistent synthetic sea surfaces from field sea states. One-dimensional and directional two-dimensional sea surface elevations are generated for each observed sea state using a spectral random-phase approach, enabling time-domain visualisation and further analysis without attempting full wave-field reconstruction.
• A methodology aligned with operational monitoring practice. The proposed approach is specifically tailored to the characteristics of real berth monitoring datasets and provides a practical bridge between raw field measurements and engineering-relevant characterisation of responses.

2. Materials and Methods

This section provides an overview of the data-processing and analysis methodology used to derive sea-state-conditioned motion response envelopes from heterogeneous berth-monitoring data. The proposed workflow is designed to accommodate incomplete overlap between motion DoF, multiple vessels, and multiple berth locations, while maintaining transparency and reproducibility.

The overall methodology consists of five main stages: data ingestion and structuring, quality control, sea-state binning, motion response envelope estimation, and synthetic sea surface generation. These stages are applied independently to each motion degree of freedom and are summarised schematically in Figure 1.

First, raw hourly records are ingested from comma-separated value files and organised into degree-of-freedom-specific tables. No assumption is made about the simultaneity of motions, vessels, or berth locations. Each record is treated as an independent observation, characterised by sea-state descriptors and corresponding motion-response statistics.

Second, a two-stage quality control procedure is applied to each degree-of-freedom dataset, as detailed in Section 2.4. This stage removes physically implausible records and statistically inconsistent response–sea-state pairs using robust regression and outlier detection. Quality control constitutes the first filtering stage of the workflow and ensures that subsequent analyses are based on physically meaningful and statistically representative data.

Third, the quality-controlled records are grouped into discrete bins defined in the space of sea-state parameters. Binning is performed with respect to significant wave height, peak wave period, and wave directional information. This step enables conditional statistical analysis without assuming stationarity of the wave field or deterministic relationships between responses.

Fourth, motion response envelopes are derived independently for each degree of freedom by computing representative and conservative statistics within each sea-state bin. These envelopes provide a compact and engineering-relevant characterisation of berth motion severity as a function of environmental conditions.

Finally, synthetic sea surface elevations corresponding to each observed sea state are generated using a spectral random-phase approach. One-dimensional time series and directionally spread two-dimensional sea surfaces are produced to facilitate time-domain interpretation, visualisation, and further response analysis. The synthetic seas are statistically consistent with the observed sea-state parameters, but do not attempt to reconstruct the historical wave field.

2.1. Dataset Description

The dataset analysed in this study is publicly available and comprises long-term field measurements of ship motions at berth collected from multiple vessels and port locations [18]. The records are reported at hourly resolution and include motion statistics together with corresponding environmental and sea-state parameters.

Motion measurements encompass combinations of surge, sway, heave, roll, pitch, and yaw, depending on vessel instrumentation and monitoring configuration. Consequently, the number of available hourly observations varies between motion degrees of freedom, as summarised in Table 3.

For each hourly record, the dataset provides key sea-state descriptors, including significant wave height $H_s$ and peak wave period $T_p$, together with wind speed and direction. Additional wave metrics are measured independently and external to the ship motions. Ship motion responses are reported as hourly mean, maximum, and standard deviation values for each available degree of freedom.

2.2. Analysis Philosophy

Given the heterogeneous nature of the dataset and the limited degree of simultaneity among motion degrees of freedom, a DoF-specific analysis strategy is adopted. Each motion component is analysed independently as a statistical response to the prevailing sea state, rather than attempting to reconstruct a fully coupled multi-DoF response. This approach maximises data utilisation while avoiding restrictive assumptions that are difficult to justify for operational berth measurements.

Hourly records are treated as independent statistical observations of the local sea state and associated vessel response. Motion statistics are conditioned directly on observed sea-state descriptors without attempting deterministic reconstruction of the wave field.

2.3. Data Loading and Pre-processing

Raw data are imported and organised using a unified processing pipeline implemented in MATLAB. For each motion DoF, the corresponding records are extracted and stored in a structured table format. Time stamps are standardised to hourly resolution, and motion statistics and sea-state parameters are associated with each hourly record.

To ensure consistency across all DoF, variable naming and formatting are harmonised during loading. Records with missing timestamps or incomplete essential fields are retained at this stage and handled subsequently during quality control, rather than being discarded a priori.

2.4. Quality Control Procedure

A two-stage quality control (QC) procedure is applied independently to each motion DoF dataset to ensure physical plausibility, internal consistency, and robustness against spurious measurements. The overall QC workflow is summarised in the methodological flowchart in Section 2, where QC serves as the first filtering stage before response envelope estimation.

2.4.1. Basic Plausibility Screening

Each hourly record is first subjected to basic plausibility checks to remove incomplete or physically unrealistic observations. Records are excluded if any of the following conditions are met:

• missing or additionally invalid time stamp;
• undefined, negative, or non-physical sea-state parameters;
• significant wave height below a minimum threshold, intended to exclude calm conditions dominated by sensor noise rather than wave-induced response;
• peak wave period outside a physically realistic range for gravity waves at berth;
• missing or undefined motion response statistic for the degree of freedom under consideration.

These criteria ensure that retained records correspond to meaningful wave-forced conditions and form a physically interpretable basis for subsequent statistical analysis.

The basic plausibility screening employed conservative thresholds to exclude noise-dominated or physically unrealistic sea states while retaining most meaningful wave-forced observations. A minimum significant wave height threshold is imposed to remove calm conditions in which measured ship motions are likely dominated by sensor noise or non-wave-related effects. Lower and upper bounds on peak wave period are applied to exclude spurious short-period fluctuations and unrealistically long periods arising from spectral estimation artefacts or non-wave low-frequency processes. The threshold values set in Table 1 are consistent with commonly adopted limits in wave measurement and statistical response analysis [4,8,17,19,20].

2.4.2. Robust Regression-Based Outlier Detection

Following basic screening, a second-stage statistical outlier detection is applied to identify inconsistent response–sea-state pairs. For each DoF, the hourly motion response standard deviation y is analysed as a function of the corresponding sea-state parameters using a robust linear regression model. The regression model is expressed by Equation (3).

$$y=\mathit{X}\mathit{\beta }+\epsilon ,$$

(3)

where $\mathit{X}$ is a matrix containing the sea-state predictors (including $H_s$, $T_p$, and wave direction), $\mathit{\beta }$ is the vector of regression coefficients, and $\epsilon$ represents the residual error vector. Wave direction is included in the regression model to account for directional variability in the vessel response during the outlier detection stage. This improves the identification of physically inconsistent response–sea-state pairs, even though direction is not used as a conditioning variable in the subsequent envelope construction.

The model is estimated using a robust M-estimation procedure with iteratively re-weighted least squares. In this framework, observations with large residuals are assigned reduced influence through adaptive weighting, thereby limiting the impact of extreme values on the fitted trend. This approach provides a stable baseline representation of the dominant response–sea-state relationship while preserving physically plausible high-response events. Robust regression is preferred over ordinary least squares because field measurements at berth may contain occasional sensor artefacts or anomalous response events that would otherwise bias the fitted relationship [21,22].

A linear model is adopted not as a physical representation of ship response, but as a local statistical baseline capturing the dominant monotonic relationship between wave severity and motion amplitude. Within the restricted ranges imposed by QC filtering and subsequent sea-state binning, a linear approximation provides a stable, interpretable reference while avoiding overfitting.

Let $\widehat{y}$ denote the fitted response obtained from the robust regression. The residuals are defined in Equation (4).

$$r=y-\widehat{y}.$$

(4)

To characterise residual dispersion in a manner insensitive to outliers, the Median Absolute Deviation (MAD) is employed as a robust scale estimator in Equation (5).

$$\mathrm{MAD}=\mathrm{median}\left(|r-\mathrm{median}(r\left)\right|\right).$$

(5)

For residuals that are approximately normally distributed, the MAD is related to the standard deviation $\sigma$ through Equation (6).

$$\sigma \approx 1.4826\times \mathrm{MAD},$$

(6)

where the factor $1.4826$ equals $1/{\Phi }^{-1}\left(0.75\right)$, with ${\Phi }^{-1}$ denoting the inverse cumulative distribution function of the standard normal distribution. This scaling ensures consistency between the MAD and the classical standard deviation estimator under the assumption of normality [21,22]. A record is classified as an outlier and excluded from further analysis if $\left|r\right|>3\times 1.4826\times \mathrm{MAD}$.

2.4.3. Threshold Selection and Sensitivity Considerations

The choice of a $3\times$ MAD threshold reflects a balance between robustness and data retention. For approximately normally distributed residuals, this criterion corresponds to a confidence level comparable to the classical $3\sigma$ rule, while remaining substantially less sensitive to extreme values than variance-based thresholds. Sensitivity testing conducted during method development indicated that lower thresholds led to excessive rejection of physically plausible high-response events, whereas higher thresholds allowed clearly inconsistent measurements to persist. The selected threshold, therefore, provides a conservative yet practical compromise for heterogeneous field datasets.

2.4.4. Role of Quality Control in the Analysis Pipeline

The QC procedure ensures that subsequent response envelope estimation is based on physically consistent and statistically representative observations. By applying identical QC logic independently to each motion DoF, the methodology avoids imposing simultaneity constraints across motions while maintaining a uniform standard of data integrity across vessels, berths, and monitoring configurations. This design choice is essential for effectively utilising long-term, heterogeneous berth-monitoring datasets.

2.5. Sea-State-Conditioned Response Envelopes

Following quality control, motion response envelopes are independently derived for each DoF using a nonparametric conditional binning approach in sea-state space. Let Y denote the motion response standard deviation for a given DoF, and let $(H_s,T_p)$ denote the significant wave height and peak wave period, respectively. The objective is to characterise the conditional behaviour of Y given the prevailing sea state.

2.5.1. Sea-State Binning

The continuous sea-state space is partitioned into rectangular bins defined by intervals in $H_s$ and $T_p$ as $H_s\in [H_{s,i}^{\mathrm{lo}},H_{s,i}^{\mathrm{hi}}],T_p\in [T_{p,j}^{\mathrm{lo}},T_{p,j}^{\mathrm{hi}}]$ where $i=1,\dots ,n_H$ and $j=1,\dots ,n_T$ denote the bin indices in wave height and period, respectively. For each bin $(i,j)$, the subset of observations is defined as in Equation (7).

$${\mathcal{Y}}_{ij}=\left\{Y_k|H_{s,k}\in [H_{s,i}^{\mathrm{lo}},H_{s,i}^{\mathrm{hi}}],T_{p,k}\in [T_{p,j}^{\mathrm{lo}},T_{p,j}^{\mathrm{hi}}]\right\}.$$

(7)

The number of observations in each bin is $N_{ij}=\mathrm{card}\left({\mathcal{Y}}_{ij}\right)$ where $\mathrm{card}(·)$ denotes set cardinality. The value $N_{ij}$ provides a measure of statistical reliability for the estimated response statistics in each bin.

2.5.2. Minimum Bin Population Requirement

To ensure statistical reliability of the estimated conditional response statistics, a minimum population threshold is imposed for each sea-state bin. Specifically, percentile-based envelope statistics are computed only if $N_{ij}\ge N_{min}$ where $N_{min}$ denotes a prescribed minimum sample size. This requirement mitigates instability in empirical quantile estimation arising from sparse data. For small sample sizes, percentile estimates are highly sensitive to individual observations and may not accurately reflect the underlying conditional distribution. By enforcing a minimum bin population, the resulting envelopes reflect statistically meaningful response behaviour while avoiding spurious amplification due to isolated or poorly sampled sea-state combinations. This threshold represents a trade-off between resolution in sea-state space and robustness of statistical inference.

2.5.3. Conditional Statistical Characterisation

For each populated bin $(i,j)$ with sufficient sample size, the empirical conditional distribution of Y given $(H_s,T_p)$ is approximated by the sample set ${\mathcal{Y}}_{ij}$. Representative and conservative statistics are then computed directly from this empirical distribution.

The conditional mean and median are given by Equations (8) and (9), respectively.

$${\overline{Y}}_{ij}={\displaystyle \frac{1}{N_{ij}}}\sum _{Y_k\in {\mathcal{Y}}_{ij}}{Y}_k,$$

(8)

$$Y_{ij}^{\left(50\right)}=Q_{0.50}\left(Y\mid H_s,T_p\right),$$

(9)

where $Q_p(·)$ denotes the empirical p-quantile operator. Upper-percentile envelopes are defined as in Equation (10) while the extreme observed value within the bin is $Y_{ij}^{\max }=max{\mathcal{Y}}_{ij}$.

$$Y_{ij}^{\left(p\right)}=Q_p\left(Y\mid H_s,T_p\right),0.90\le p\le 0.95$$

(10)

2.5.4. Definition of the Response Envelope

The sea-state-conditioned response envelope is defined as the mapping.

$${\mathcal{E}}_p(H_s,T_p)=Q_p\left(Y\mid H_s,T_p\right),$$

(11)

where p represents the selected percentile level, the median envelope ($p=0.50$) describes the typical response behaviour. In contrast, upper-percentile envelopes ($0.90\le p\le 0.95$) represent conservative response levels relevant to berth operability assessment. The upper-percentile envelopes are defined as such, providing a statistically stable representation of high-response behaviour while avoiding excessive sensitivity to individual extreme observations within sparsely populated sea-state bins.

This formulation constitutes a non-parametric estimation of the conditional distribution $f_{Y|H_s,T_p}(y\mid h,t)$ without imposing assumptions regarding linearity, normality, or specific response transfer functions. The envelope, therefore, reflects empirically observed berth response behaviour across heterogeneous vessels and locations while preserving statistical robustness.

Directional values are retained across the full angular range (0 – 360°), and no directional filtering is applied during the response-envelope estimation. Motion statistics are therefore conditioned only on the sea-state parameters $H_s$ and $T_p$. Wave direction is analysed separately to assess directional trends in the response, but is not treated as a conditioning variable in the envelope construction. This approach preserves the full dataset while avoiding excessive subdivision of the observations into sparsely populated directional bins in the multi-dimensional sea-state space.

2.6. Synthetic Sea Surface Generation

For each quality-controlled hourly record, synthetic sea surface elevations are generated to provide statistically consistent time-domain realisations of the observed sea states. The objective is not to reconstruct the historical wave field, but to produce representative surface elevation fields consistent with the measured sea-state descriptors for visualisation, interpretation, and further response analysis.

2.6.1. Numerical Synthesis Settings

The numerical parameters adopted for sea surface synthesis are summarised in Table 2.

For the one-dimensional formulation, a time step of $0.5\mathrm{s}$ and a duration of one hour are selected to ensure adequate temporal resolution while maintaining consistency with the hourly response statistics. The spectrum is discretised using 512 frequency components, providing sufficient resolution of the prescribed peak period and spectral bandwidth without introducing artificial periodicity. A JONSWAP spectrum [23] with peak enhancement factor $\gamma$ value of 3.3 is adopted, consistent with widely accepted engineering practice for wind-generated seas [1,8,24]. The parameter $\gamma$ governs the sharpness of the spectral peak. When this parameter value is set to one, it reduces the spectrum to the Pierson–Moskowitz form [25], whereas larger values represent developing wind seas with narrower spectral concentration. Sensitivity analysis confirmed that moderate variations in $\gamma$ primarily affect the spectral bandwidth and the visual appearance of the realisations, while preserving the prescribed significant wave height and peak period.

For the two-dimensional directional synthesis, the objective is illustrative rather than statistical reconstruction. The duration is therefore limited to 10 minutes, with a temporal resolution of $1.0\mathrm{s}$ to reduce computational cost while adequately resolving dominant wave periods. Since crest geometry and spatial coherence are governed primarily by spectral discretisation and spatial resolution, this shorter duration does not compromise physical realism. The resulting 600-frame realisations provide representative visualisations of directional spreading consistent with the prescribed sea-state parameters.

Wave direction is obtained directly from the dataset and is used to define the dominant propagation direction in the directional spectrum for the synthetic sea surface generation. No attempt is made to infer wave direction from ship motions, since reliable directional information is already available from the environmental measurements. The directional parameter is therefore treated as an observed sea-state descriptor rather than an estimated quantity. In harbour environments influenced by reflection, diffraction, and local geometry, the measured wave direction should be interpreted as an approximate indicator of dominant forcing orientation rather than a single deterministic incident direction. The adopted approach maintains consistency with the statistical nature of the response-envelope framework while avoiding unnecessary inverse estimation of wave direction.

The selected synthesis parameters follow established spectral wave-modelling practice and are consistent with recommendations in Faltinsen [1], Holthuijsen [8], and Goda [24]. Variations in discretisation and spreading parameters primarily affect the visual texture of the synthetic sea surfaces and do not alter the enforced low-order statistical consistency with the prescribed sea-state descriptors.

2.6.2. Spectral Formulation and Normalisation

The target frequency spectrum is modelled using a Pierson–Moskowitz (PM) base shape with JONSWAP peak enhancement [23,25]. Let ${\omega }_p=2\pi /T_p$ denote the peak angular frequency. The defined unscaled spectral shape is given by Equation (12).

$$\tilde{S}\left(\omega \right)={\omega }^{-5}exp\left[-{\displaystyle \frac{5}{4}}{\left({\displaystyle \frac{{\omega }_p}{\omega }}\right)}^4\right]{\gamma }^{exp\left[-{\displaystyle \frac{{(\omega -{\omega }_p)}^2}{2{\kappa }^2{\omega }_p^2}}\right]}$$

(12)

where $\kappa$ is a parameter that controls the peak width.

Rather than relying on closed-form amplitude constants, the spectrum is numerically normalised so that its discrete zeroth moment matches the prescribed significant wave height. The zeroth spectral moment of the wave spectrum is related to the significant wave height, as shown in Equation (17), a standard result in spectral wave theory [4,8]. Using a discrete frequency grid ${\left\{{\omega }_n\right\}}_{n=1}^{N_{\omega }}$ with spacing $\Delta \omega$, the raw moment is calculated by Equation (13).

$$m_0^{\mathrm{raw}}=\sum _{n=1}^{N_{\omega }}\tilde{S}\left({\omega }_n\right)\Delta \omega$$

(13)

The scaling factor is then defined by Equation (14).

$$A={\displaystyle \frac{m_0}{m_0^{\mathrm{raw}}}}$$

(14)

so that the final spectrum used in synthesis is

$$S\left({\omega }_n\right)=A\tilde{S}\left({\omega }_n\right).$$

(15)

This normalisation preserves the spectral shape governed by $(T_p,\gamma )$ while enforcing consistency with $H_s$ through discrete integration. The same scaled spectrum $S\left(\omega \right)$ is used in both the one-dimensional and directional two-dimensional sea surface synthesis.

2.6.3. One-Dimensional Sea Surface Representation

The one-dimensional sea surface elevation $\eta \left(t\right)$ is synthesised using a linear random-phase superposition of harmonic wave components. The surface elevation is expressed as Equation (16). This linear random-phase representation follows the classical spectral description of irregular waves introduced by Longuet-Higgins [19], in which the sea surface elevation is modelled as a Gaussian random process obtained by superposing harmonic components with statistically independent phases.

$$\eta \left(t\right)=\sum _{n=1}^N\sqrt{2S\left({\omega }_n\right)\Delta \omega }cos\left({\omega }_{n}t+{\varphi }_n\right),$$

(16)

where $S\left(\omega \right)$ is the wave energy spectrum, ${\omega }_n$ denotes the angular frequency of the n-th spectral component, $\Delta \omega$ is the frequency increment, and ${\varphi }_n$ are independent random phases uniformly distributed on the range $[0,2\pi ]$.

Equation (16) represents a specific case of harmonic superposition in which the component amplitudes are determined directly from the prescribed wave spectrum. In contrast to a general Fourier representation with arbitrary amplitudes, the choice $a_n=\sqrt{2S\left({\omega }_n\right)\Delta \omega }$ ensures that the variance of the synthesised surface elevation is consistent with the zeroth spectral moment, thereby preserving the significant wave height and spectral energy distribution.

The frequency range $[{\omega }_{min},{\omega }_{max}]$ and the number of spectral components N are selected to ensure adequate resolution of the energy-containing portion of the spectrum and numerical consistency of the discrete spectral integration. The upper frequency limit is constrained by the Nyquist criterion associated with the chosen time step $\Delta t$, ensuring that ${\omega }_{max}\le \pi /\Delta t$ to avoid temporal aliasing. The lower frequency bound and spectral discretisation are chosen to capture the dominant energy near the spectral peak and to ensure an accurate approximation of the zeroth spectral moment $m_0$. The spectral density $S\left(\omega \right)$ is parameterised using the observed significant wave height $H_s$ and peak wave period $T_p$, illustrated in Equations (17) and (18). The resulting realisation $\eta \left(t\right)$ represents a stationary, Gaussian sea surface consistent with the target spectral parameters over the synthesis interval.

$$H_s=4\sqrt{m_0}$$

(17)

where $m_0$ is the zeroth spectral moment described in Equation (18).

$$\begin{array}{cc}\hfill m_0& ={\int }_0^{\infty }S\left(\omega \right)\mathrm{d}\omega \hfill \\ \hfill & ={\displaystyle \frac{H_s^2}{16}}\hfill \end{array}$$

(18)

2.6.4. Directional Two-Dimensional Sea Surface Representation

To account for wave directionality, a two-dimensional directional sea surface $\eta (x,y,t)$ is generated by extending the spectral formulation to include directional spreading. The surface elevation is expressed by Equation (19).

$$\eta (x,y,t)=\sum _{n=1}^{N_{\omega }}\sum _{m=1}^{N_{\theta }}\sqrt{2S({\omega }_n,{\theta }_m)\Delta \omega \Delta \theta }cos\left({\mathbf{k}}_{n,m}·\mathbf{x}-{\omega }_{n}t+{\varphi }_{n,m}\right),$$

(19)

where $\mathbf{x}=(x,y)$ is the horizontal position vector, ${\mathbf{k}}_{n,m}$ is the wave vector corresponding to frequency ${\omega }_n$ and direction ${\theta }_m$, and ${\varphi }_{n,m}$ are independent random phases. The frequency spectrum $S\left(\omega \right)$ defined in Section 2.6.2 is extended to a directional form through factorisation, in Equation (20).

$$S(\omega ,\theta )=S\left(\omega \right)D\left(\theta \right),$$

(20)

where $D\left(\theta \right)$ represents the directional spreading function centred around the dominant wave direction ${\theta }_0$. In this study, the observed wave direction in the dataset is used as the central direction. A cosine-power spreading model is adopted to represent directional energy distribution, which is widely used in spectral wave modelling [8,24]. The spreading function is defined in Equation (21).

$$D\left(\theta \right)=\left\{\begin{array}{cc}C{cos}^{2s}(\theta -{\theta }_0),\hfill & |\theta -{\theta }_0|\le {\displaystyle \frac{\pi }{2}},\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(21)

where ${\theta }_0$ denotes the mean wave direction and s is the directional spreading parameter controlling the concentration of wave energy around ${\theta }_0$. Larger values of s correspond to narrower directional spreading, whereas smaller values produce broader directional distributions. The normalisation constant ensures that ${\int }_{-\pi }^{\pi }D\left(\theta \right)d\theta =1$ .

The wave vector magnitude is obtained from the linear dispersion relation of surface gravity waves, given by [1,8], as in Equation (22).

$${\omega }^2=gktanh\left(kh\right),$$

(22)

where g is the gravitational acceleration and h is the local water depth. In the absence of detailed bathymetric information, deep-water conditions are assumed as a practical approximation. Under deep-water conditions ($kh\gg 1$), the hyperbolic tangent term approaches unity, i.e., $tanh\left(kh\right)\approx 1$, and the dispersion relation reduces to ${\omega }^2\approx gk$, from which the wave number is obtained as in Equation (23).

$$k={\displaystyle \frac{{\omega }^2}{g}}$$

(23)

This approximation is valid when the water depth is large relative to the wavelength (approximately $h⩾\lambda /2$), such that finite-depth effects on wave propagation are negligible.

2.7. Verification and Sensitivity of Synthetic Sea-State Reconstruction

Because the synthetic sea surfaces are used to provide time-domain realisations that are statistically consistent with the observed sea-state descriptors (rather than to reconstruct the historical wave field), verification focuses on matching low-order spectral statistics implied by $(H_s,T_p)$ and on demonstrating robustness of downstream operability metrics to synthesis settings.

2.7.1. Statistical Consistency Checks

For each generated sea realisation, the reconstructed elevation time series $\eta \left(t\right)$ is verified against the target variance implied by the significant wave height. For a linear Gaussian sea state, the relationship between the significant wave height and the sample variance is depicted in Equations (17) and (18).

From the synthetic realisation, the sample variance is computed as shown in Equation (24).

$${\widehat{\sigma }}_{\eta }^2={\displaystyle \frac{1}{N_t-1}}\sum _{k=1}^{N_t}{\left(\eta \left(t_k\right)-\overline{\eta }\right)}^2,$$

(24)

The agreement between the synthetic realisation and the sample variances is quantified using a relative error calculated using Equation (25).

$${ϵ}_{\sigma }={\displaystyle \frac{\left|16{\widehat{\sigma }}_{\eta }^2-{H_s}^2\right|}{{H_s}^2}}$$

(25)

Additionally, the sample spectrum $\widehat{S}\left(\omega \right)$ estimated from $\eta \left(t\right)$ is compared to the target $S\left(\omega \right)$ to confirm that energy is concentrated near the prescribed peak frequency ${\omega }_p=2\pi /T_p$. A practical diagnostic is the location of the dominant spectral peak ${\widehat{\omega }}_p$, and the peak-period mismatch is quantified in Equation (26).

$${ϵ}_p={\displaystyle \frac{|{\widehat{T}}_p-T_p|}{T_p}}$$

(26)

These checks ensure that each synthetic realisation is consistent with the sea-state descriptors used to condition the motion response envelopes.

2.7.2. Ablation Study: Sensitivity to Synthesis Parameters

An ablation study is performed to verify that the statistical envelope results and derived operability metrics are not unduly sensitive to numerical synthesis settings. In addition to envelope-based motion metrics, the present ablation explicitly evaluates the low-order consistency of the synthesised sea states with their target parameters. Let $\mathcal{M}$ denote any scalar performance or consistency metric extracted from the synthesised realisations.

In particular, two primary consistency metrics are considered:

1.

Relative error in significant wave height. Let $H_s$ denote the target significant wave height for a given sea-state bin, and let ${\widehat{H}}_s$ denote the value re-estimated from the synthesised time series $\eta \left(t\right)$ using the spectral moment method. The relative error is defined by Equation (27).

$${\mathcal{M}}_{H_s}={\displaystyle \frac{\left|{\widehat{H}}_s-H_s\right|}{max(H_s,\epsilon )}}$$

(27)

where $\epsilon$ is a small numerical floor threshold (e.g., ${10}^{-6}$) introduced to prevent division by zero in near-calm sea states.

2.

Relative error in peak period. Let $T_p$ denote the target peak period of the JONSWAP spectrum and ${\widehat{T}}_p$ the peak period estimated from the reconstructed spectrum of the synthesised signal. The relative error is defined in Equation (28).

$${\mathcal{M}}_{T_p}={\displaystyle \frac{\left|{\widehat{T}}_p-T_p\right|}{max(T_p,\epsilon )}}$$

(28)

These metrics quantify the degree to which the synthesised realisations preserve the prescribed first- and second-order spectral characteristics of the sea state. They are referred to here as operability consistency metrics because they do not measure motion-limit exceedance directly, but instead verify that the underlying stochastic forcing statistics are faithfully reproduced; reliable envelope-based operability conclusions depend on this low-order statistical consistency.

For a baseline synthesis configuration ${\theta }_0$ and a perturbed configuration $\theta$, sensitivity is quantified as in Equation (29).

$$\Delta \mathcal{M}\left(\theta \right)=\mathcal{M}\left(\theta \right)-\mathcal{M}\left({\theta }_0\right),$$

(29)

evaluated across representative sea-state bins $(H_s,T_p)$.

The ablation varies (one-at-a-time) the key parameters controlling discretisation and directional spreading, including the number of frequency components $N_{\omega }$, time step $\Delta t$, synthesis duration T, the JONSWAP peak enhancement factor $\gamma$, and for 2D fields the number of directional components $N_{\theta }$ and spreading parameter s. The intent is to show that these parameters influence the visual texture and sampling variability of the realisations, while preserving the enforced low-order consistency with $(H_s,T_p)$ and leaving envelope-based conclusions unchanged within the expected Monte Carlo variability.

2.7.3. Why this Validation is Appropriate for the Present Framework

The proposed framework does not require deterministic reconstruction of harbour wave directionality or phase history. Instead, it requires that synthetic seas be statistically consistent with the descriptors used for conditioning the observed motion envelopes. The verification targets ${\sigma }_{\eta }^2$ (via $H_s$), peak location (via $T_p$), and spectral energy concentration, which are the quantities directly implied by the measured sea-state parameters and used to interpret envelope-based response levels in the time domain. This design provides a transparent check on internal consistency while remaining aligned with the study objective of field-driven response characterisation.

3. Results

This section presents the results of applying the proposed sea-state-conditioned analysis framework to the berth motion dataset. Results are reported separately for each motion degree of freedom. They are structured to reflect the sequence of analysis steps: quality-controlled data availability, sea-state-conditioned response envelopes, directional effects, and summary observations.

3.1. Quality-Controlled Data Availability

Application of the two-stage quality control procedure described in Section 2.4 resulted in substantial filtering of the raw hourly records. The number of retained observations varies markedly across motion DoF, see Table 3, reflecting differences in vessel instrumentation, berth locations, and monitoring durations.

Table 3. Number of hourly records available for each motion degree of freedom before and after quality control.

DoF	Initial Records	After Plausibility Screening	After Outlier Removal
Surge	361	361	344
Sway	1242	1242	1188
Heave	359	359	357
Roll	1302	1302	1240
Pitch	1302	1302	1150
Yaw	1225	1225	1091

Table 3. Number of hourly records available for each motion degree of freedom before and after quality control.

DoF	Initial Records	After Plausibility Screening	After Outlier Removal
Surge	361	361	344
Sway	1242	1242	1188
Heave	359	359	357
Roll	1302	1302	1240
Pitch	1302	1302	1150
Yaw	1225	1225	1091

Despite limited temporal overlap among motions, each DoF dataset retained sufficient records after quality control to enable statistically meaningful conditional analyses. By treating each motion independently, the proposed framework avoids the significant data loss that would arise from enforcing simultaneity across the DoF.

Table 3 summarises the record availability for each motion degree of freedom before and after quality control. Figure 2 presents the corresponding response–sea-state scatter plots. Retained observations are shown in blue, while points removed during plausibility screening (grey) and robust outlier detection (red) are also indicated.

3.2. Response–Sea-State Relationships

Figure 2 presents the quality-controlled response–sea-state scatter plots for all motion degrees of freedom. The vertical axis represents hourly motion variability, defined as the standard deviation of the motion response within each hourly record.

A consistent positive dependence between significant wave height $H_s$ and motion variability is observed across all motions, reflecting increased excitation under more energetic sea states. While this overall monotonic trend is evident across all degrees of freedom, the degree of dispersion varies. Vertical motions exhibit comparatively compact response bands, whereas horizontal and rotational motions display broader scatter, indicating greater sensitivity to berth geometry, mooring configuration, and directional forcing.

The quality-control procedure preserves the dominant response–sea-state structure while removing isolated high-residual or physically inconsistent observations. In particular, robust regression–based outlier detection reduces excessive dispersion at moderate $H_s$ levels without distorting the underlying trend. This improves the dataset’s statistical coherence and provides a stable basis for subsequent envelope estimation.

3.3. Sea-State-Conditioned Motion Response Envelopes

Sea-state-conditioned motion response envelopes are derived by binning the quality-controlled records in the $(H_s,T_p)$ space. For each bin, representative and conservative statistics of the motion standard deviation are computed, yielding empirical response envelopes conditioned on the observed sea state. The resulting envelopes are summarised in Figure 3 and provide a compact statistical characterisation of berth motion behaviour without requiring response amplitude operators or detailed vessel and mooring models.

Across all degrees of freedom, the median envelopes exhibit a systematic increase with increasing significant wave height, indicating that wave energy is the dominant driver of berth motions.

3.3.1. Quantitative Characterisation of Motion Response

To provide a quantitative comparison between motion components, additional summary metrics are derived from the envelope curves. These metrics characterise (i) the sensitivity of each motion to increasing wave height, (ii) the variability between typical and extreme responses, and (iii) the normalised motion response relative to wave height.

First, a motion sensitivity coefficient is estimated from the slope of the median response envelope with respect to significant wave height. This coefficient quantifies how strongly motion variability increases with increasing wave energy as depicted by Equation (30).

$$S_m={\displaystyle \frac{\Delta {\sigma }_{\mathrm{motion}}}{\Delta H_s}}$$

(30)

Second, the envelope spread is quantified as the difference between the upper-percentile (P90) envelope and the median envelope at a given wave height is given by Equation (31).

$$\Delta {\sigma }_{\mathrm{motion}}\left(H_s\right)={\sigma }_{\mathrm{motion},P90}\left(H_s\right)-{\sigma }_{\mathrm{motion},50}\left(H_s\right)$$

(31)

where ${\sigma }_{\mathrm{motion},P90}$ and ${\sigma }_{\mathrm{motion},50}$ denote the median and 90th percentile response levels within each sea-state bin.

Finally, a motion coefficient is introduced in Equation (32) to normalise the response by wave height.

$$C_m={\displaystyle \frac{{\sigma }_{\mathrm{motion}}}{H_s}}$$

(32)

This coefficient provides a statistical analogue of a response amplitude operator and enables direct comparison between motions with different physical units and magnitudes.

Table 4 summarises the resulting mean quantitative metrics for all six degrees of freedom.

The quantitative metrics reveal clear differences in the response characteristics of the various motion components. Sensitivity values indicate that roll and surge exhibit the strongest dependence on wave height, with sensitivity coefficients of approximately 0.22 and 0.21, respectively. Heave shows moderate sensitivity, whereas sway, yaw, and pitch respond more weakly to increasing wave height.

Envelope spread indicates the variability between typical and extreme response levels. Roll displays the largest spread (0.44), followed by yaw (0.34), suggesting that these rotational motions are particularly susceptible to intermittent amplification. Surge and sway exhibit moderate spreads, while pitch and heave show comparatively smaller values.

The motion coefficient further highlights the relative magnitude of motions when normalised by wave height. Roll exhibits the largest coefficient ($C_m\approx 0.39$), indicating strong amplification relative to wave height and reflecting the well-known sensitivity of roll to resonant excitation [1,4,5,6]. Surge displays the largest coefficient among the translational motions, consistent with the influence of low-frequency horizontal forcing.

Pitch shows both the lowest sensitivity and the smallest envelope spread, suggesting comparatively stable behaviour across the analysed sea states. Yaw, in contrast, exhibits moderate sensitivity but a large envelope spread, indicating that extreme yaw responses occur intermittently under particular combinations of wave period and directional forcing.

3.3.2. Quadratic Envelope Fits and Curvature of Motion Response

To quantify the curvature in the envelope plots, quadratic functions were fitted separately to the median and P90 response envelopes for each degree of freedom as functions of significant wave height $H_s$, according to Equation (33).

$${\sigma }_{\mathrm{mov}}\left(H_s\right)=a_0+a_1{H}_s+a_2{H}_s^2$$

(33)

The quadratic coefficient $a_2$ provides a direct measure of curvature, while the goodness of fit is summarised by the coefficient of determination $R^2$. The fitted curves are shown in Figure 3, and the corresponding coefficients are summarised in Table 5. Only the quadratic curvature coefficient and the corresponding goodness-of-fit are reported here, as these quantities are the most directly relevant for interpreting superlinear growth of the response envelopes. The full fitted coefficients are provided in Appendix B.

The fitted results confirm that the motion envelopes are not purely linear in $H_s$, although the strength and regularity of the curvature differ between motion components. For the median envelopes, roll exhibits the strongest positive quadratic coefficient ($a_{2,\mathrm{med}}=0.0376$, $R^2=0.994$), followed by yaw ($a_{2,\mathrm{med}}=0.0223$, $R^2=0.984$) and surge ($a_{2,\mathrm{med}}=0.0139$, $R^2=0.887$). This indicates clear superlinear growth in the typical response of these motions with increasing wave height. Pitch and sway also show positive but weaker median curvature. In contrast, heave is the only motion with a slightly negative quadratic coefficient in the median fit ($a_{2,\mathrm{med}}=-0.0101$, $R^2=0.976$), consistent with its comparatively smooth and nearly linear behaviour over the analysed $H_s$ range.

For the P90 envelopes, the strongest positive curvature is observed for heave ($a_{2,\mathrm{P}90}=0.0257$, $R^2=0.914$), roll ($a_{2,\mathrm{P}90}=0.0212$, $R^2=0.841$), yaw ($a_{2,\mathrm{P}90}=0.0168$, $R^2=0.746$), and sway ($a_{2,\mathrm{P}90}=0.0145$, $R^2=0.991$). These results indicate that the conservative response levels for these motions increase more rapidly than linearly with wave height. In contrast, the P90 surge fit has a negative quadratic coefficient ($a_{2,\mathrm{P}90}=-0.0338$, $R^2=0.759$), which reflects the abrupt intermediate-$H_s$ jump and subsequent partial flattening seen in the plotted envelope rather than a smooth monotonic superlinear trend. The lower $R^2$ values for surge and yaw P90 also suggest that the upper-tail behaviour of these motions is more irregular and less well represented by a simple quadratic model.

A useful distinction emerges between typical and extreme response curvature. Roll, surge, and yaw show clear curvature in the median envelope, indicating that their typical responses strengthen nonlinearly with increasing sea severity. By contrast, heave and sway display stronger curvature in the P90 envelope than in the median envelope, suggesting that their extreme responses amplify more rapidly than their central tendency. Pitch shows relatively weak curvature in both envelopes, which is consistent with the previously observed modest sensitivity and limited spread of this motion component.

Overall, the quadratic fits show that envelope curvature is motion-dependent and that nonlinear growth with $H_s$ is more pronounced for some response modes than others. Roll emerges as the clearest example of systematic superlinear growth across both typical and conservative response levels, while yaw and surge exhibit strong curvature alongside greater upper-tail irregularity. Heave remains comparatively regular in the median response, but its upper-percentile envelope still shows appreciable curvature. These results reinforce the value of separating median and upper-percentile envelopes, since the nonlinear growth of extreme response levels is not always captured by the median trend alone.

3.4. Directional Dependence of Motion Response

The influence of wave direction on berth motions is examined by grouping the quality-controlled records into discrete directional sectors based on the observed wave propagation direction. For each sector, the distribution of the motion standard deviation $\mathrm{mov}\text{\_}\mathrm{sig}$ is summarised using boxplots.

Figure 4 presents the sector-wise response distributions for all six degrees of freedom. Directional variability refers to systematic changes in the conditional distribution of motion response across wave-direction sectors, manifested as sector-to-sector differences in median response level, spread, or upper-percentile behaviour.

The results indicate that several motion components exhibit discernible directional structure. In particular, roll and yaw display pronounced sector-to-sector variability. For roll, the median response increases from approximately $0.6$ in the 0–${30}^{\circ }$ sector to roughly $1.1$–$1.4$ in the 300–${360}^{\circ }$ sectors, with upper values approaching $2.5$–$3.0$. A similar but less pronounced pattern is observed for yaw, where median responses increase from about $0.12$–$0.15$ in the 0–${30}^{\circ }$ sector to approximately $0.30$–$0.35$ in the 300–${360}^{\circ }$ sectors, with several upper values exceeding $1.0$. These sector-to-sector differences suggest strong sensitivity to the orientation of incident waves relative to the vessel and berth geometry.

Sway and pitch show more moderate directional dependence. For sway, the median response increases from roughly $0.25$ in the 270–${300}^{\circ }$ sector to about $0.65$–$0.7$ in the 0–${30}^{\circ }$ sector, while pitch medians increase from approximately $0.15$ in the 0–${30}^{\circ }$ sector to around $0.28$–$0.30$ in the 300–${360}^{\circ }$ sectors. In contrast, the vertical motion component (heave) exhibits comparatively weaker directional sensitivity, with median values increasing only from approximately $0.20$ in the 0–${30}^{\circ }$ sector to about $0.35$–$0.40$ in the 300–${360}^{\circ }$ sectors.

These observations are consistent with the physical mechanisms governing berth dynamics. Rotational and horizontal motions are strongly influenced by the alignment of incoming waves relative to the vessel heading and mooring configuration. In contrast, vertical motion is primarily controlled by wave amplitude and is less sensitive to directional variations.

Note that the markers shown outside the whiskers in Figure 4 correspond to standard boxplot outliers based on the interquartile range criterion and do not indicate observations removed during quality control. All plotted data satisfy the quality-control criteria and represent valid, though extreme, response events within the retained dataset.

3.5. Synthetic Sea Surface Illustrations

To support time-domain interpretation of the statistically derived response envelopes, synthetic sea surface elevations are generated for representative quality-controlled sea states. Figure 5 presents an illustrative example consisting of a one-dimensional sea surface elevation time series and corresponding two-dimensional spatial snapshots extracted at a selected time instant.

Figure 5(a) shows the one-dimensional sea surface elevation $\eta \left(t\right)$ generated for a representative sea state. The time series exhibits irregular wave behaviour consistent with a random-phase spectral representation, with crest and trough amplitudes compatible with the prescribed significant wave height. Although no attempt is made to reproduce the historical wave record, the temporal variability provides a physically plausible reference for interpreting motion response statistics in the time domain.

Figure 5(b) and (c) present two-dimensional snapshots of the sea surface elevation $\eta (x,y,t)$ at the same time instant. The contour plot highlights the spatial distribution of wave elevation across the computational grid. At the same time, the three-dimensional surface view provides an intuitive visualisation of wave crest geometry and relative amplitude. The spatial patterns reveal coherent wave structures aligned with the prescribed wave direction, together with spatial variability arising from the superposition of multiple frequency and directional components.

The spatial extent and grid resolution of the two-dimensional domain are selected to balance visual clarity with computational efficiency. They are not intended to represent the full spatial complexity of the berth environment. Instead, these snapshots serve to illustrate the statistically consistent wave fields associated with the sea states used to derive the response envelopes.

By linking envelope-based motion statistics to representative synthetic sea surface realisations, the results provide an intuitive connection between abstract statistical characterisation and physically interpretable wave fields. This connection facilitates qualitative assessment of motion exceedance behaviour. It supports the use of envelope-derived response levels in berth operability evaluation without requiring deterministic reconstruction of site-specific wave conditions.

3.6. Summary of Observed Response Characteristics

The results demonstrate that the proposed framework successfully extracts physically meaningful motion response characteristics from heterogeneous field data. Sea-state-conditioned response envelopes provide a compact and interpretable representation of berth motion severity across a wide range of environmental conditions. The independent treatment of each degree of freedom maximises data utilisation while maintaining a consistent quality standard across vessels and berth locations.

3.7. Mini-Ablation Study of Sea Synthesis Parameters

To evaluate the numerical robustness and statistical consistency of the synthetic sea surface generation procedure, a targeted mini-ablation study is conducted. Key synthesis parameters are systematically varied in both the 1D and 2D formulations. The reconstructed significant wave height $H_s$ and peak wave period $T_p$ are compared with their prescribed values, using median relative error as the performance metric.

Table 6 summarises the median relative errors obtained for all tested configurations, while Figure 6 illustrates the full error distributions.

3.7.1. One-Dimensional Sea Synthesis

For the 1D random-phase synthesis, the baseline configuration reproduced prescribed sea-state parameters with sub-per cent accuracy, yielding median relative errors of $0.28\%$ for $H_s$ and $0.82\%$ for $T_p$. This confirms strong statistical convergence of the linear spectral superposition.

Simulation duration had the strongest influence on $H_s$. Reducing the duration to 1800 s increased the median $H_s$ error to $0.94\%$, whereas extending it to 7200 s reduced the error to $0.23\%$. This behaviour reflects the stochastic convergence properties of variance-based statistics, since $H_s=4\sqrt{m_0}$ and $m_0={\int }_0^{\infty }S\left(\omega \right)\mathrm{d}\omega$. Longer time series improve the estimation of $m_0$ and, consequently, of $H_s$.

The peak enhancement factor $\gamma$ primarily affected $T_p$. Reducing $\gamma$ to 1.0 increased the median $T_p$ error to $2.66\%$, consistent with spectral peak broadening. Increasing $\gamma$ to 5.0 reduced the error to $0.65\%$.

Frequency resolution, defined here as the number of spectral frequency components $N_{\omega }$ used to discretise the continuous wave spectrum, produced moderate variation in $T_p$ reconstruction accuracy. Median relative errors ranged from $1.06\%$ for $N_{\omega }=256$ to $1.59\%$ for $N_{\omega }=1024$, while the reconstructed $H_s$ remained below approximately $0.45\%$ in all cases.

These results indicate that peak-period estimation is moderately sensitive to spectral discretisation density. In contrast, the zeroth spectral moment (and hence $H_s$) is comparatively robust to changes in $N_{\omega }$ within the tested range.

Variation of the temporal sampling interval $\Delta t$ had negligible influence on either $H_s$ or $T_p$, confirming that the adopted time resolution satisfies the Nyquist criterion for the energy-containing portion of the spectrum.

Overall, the 1-D formulation demonstrates high numerical stability and limited sensitivity to moderate variations in parameters.

3.7.2. Two-Dimensional Sea Synthesis

For the directional 2D synthesis, errors are evaluated at a representative spatial grid point. The baseline configuration yielded median relative errors of $10.2\%$ for $H_s$ and $4.22\%$ for $T_p$.

These deviations are expected, as directional superposition produces spatial interference patterns, so the local variance at a single point does not necessarily equal the spatially integrated spectral variance.

Frequency resolution is identified as the dominant control parameter. Reducing the number of frequency components to 64 increased the median $H_s$ error to $20.9\%$, whereas increasing it to 256 reduced the error to $5.08\%$. Corresponding $T_p$ errors ranged from $3.01\%$ to $5.97\%$.

In contrast, directional discretisation had limited influence on $H_s$ (10.3–10.7%), and variation of the spreading parameter produced similarly modest changes (10.5–10.8%). These findings indicate that directional resolution primarily affects spatial texture rather than bulk spectral variance.

The synthetic sea surfaces are generated using random-phase superposition. Consequently, finite-duration realisations exhibit sampling variability in their estimated variance and peak period, even when the underlying target spectrum is identical.

Sensitivity testing using multiple random seeds confirmed that the observed variability in reconstructed $H_s$ for the 2D case is primarily attributable to stochastic phase realisation effects rather than numerical instability of the spectral formulation. The variance fluctuations decrease with increasing simulation duration and frequency resolution, consistent with Monte-Carlo convergence behaviour of random-phase spectral synthesis. Representative multi-seed results are provided in Appendix A.

3.7.3. Implications for Parameter Selection

The ablation study confirms that the adopted baseline synthesis parameters provide a stable compromise between computational efficiency and statistical fidelity. The 1D formulation reproduces prescribed sea-state descriptors with sub-per cent accuracy, while the 2D formulation produces physically realistic directional fields with acceptable local variability. No variation in the tested parameters introduced systematic bias or instability. These results validate the numerical robustness of the sea synthesis framework and support its use for time-domain visualisation and envelope-based interpretation.

4. Discussion

4.1. Comparison with Existing Berth Motion Assessment Approaches

Traditional assessment of ship motions at berth has relied heavily on frequency-domain RAOs derived from linear potential-flow theory and idealised mooring models. While such approaches provide valuable insight into resonance mechanisms, their application in operational port environments is often constrained by limited knowledge of vessel particulars, mooring properties, and local wave transformation processes.

The results presented in this study demonstrate that empirically derived, sea-state-conditioned response envelopes capture consistent and physically interpretable motion behaviour across multiple vessels and berth locations without recourse to detailed deterministic modelling. In contrast to RAO-based methods, the proposed framework accommodates heterogeneous datasets and incomplete overlap between motion degrees of freedom, which are typical of long-term operational monitoring campaigns.

4.2. Interpretation of Envelope Variability

The response–sea-state relationships indicate that motion variability generally increases with significant wave height $H_s$, confirming that wave energy is the primary driver of berth motions. However, the degree of variability differs between motion components, reflecting the combined influence of hydrodynamic forcing, vessel dynamics, and mooring system behaviour.

Quantitative envelope metrics provide additional insight into these differences. Sensitivity values indicate that roll and surge exhibit the strongest dependence on wave height, with sensitivity coefficients of approximately 0.22 and 0.21, respectively. Heave shows moderate sensitivity despite being directly forced by vertical wave motion. This behaviour likely reflects the influence of berth constraints and mooring stiffness, which limit vertical motion and reduce the rate at which heave variability increases with wave height. In contrast, horizontal and rotational motions are more strongly affected by low-frequency excitation and mooring dynamics, resulting in higher sensitivity values.

Clear differences also emerge between translational and rotational motions in terms of the variability between typical and extreme response levels. Rotational motions, particularly roll and yaw, exhibit the largest variability between median and upper-percentile responses. This behaviour is consistent with the sensitivity of rotational modes to resonant excitation, vessel stability characteristics, and the interaction between wave forcing and mooring system dynamics.

Roll amplification is particularly pronounced, reflecting the well-known susceptibility of roll motion to resonance effects and nonlinear restoring forces associated with vessel stability. In contrast, vertical motion (heave) exhibits comparatively regular behaviour. Because heave is primarily governed by direct wave excitation, its response tends to scale more smoothly with wave height and is less influenced by low-frequency drift forces or mooring constraints.

Horizontal surge and sway motions occupy an intermediate regime. These responses are influenced by both wave excitation and low-frequency drift forces, which can introduce intermittent amplification under certain sea-state conditions. Such behaviour is consistent with the combined effects of long-period wave components, mooring stiffness, and berth geometry.

From an operational perspective, the distinction between typical and extreme response behaviour is particularly important. While median envelopes represent the dominant motion levels experienced at berth, upper-percentile responses capture intermittent amplification events that frequently govern berth operability limits and the safety of the mooring system.

4.3. Implications for Berth Operability and Monitoring Practice

The combination of response envelopes and statistically consistent synthetic sea surface realisations provides a practical framework for translating field measurements into operability-relevant metrics. Rather than requiring deterministic simulations for specific vessel configurations, the proposed approach enables screening-level assessment of motion limits, exceedance likelihood, and relative berth performance under varying sea states.

This framework is particularly suited to ports where long-term monitoring data exist, but detailed vessel and mooring information is unavailable or impractical to obtain. By focusing on statistically robust field-observed behaviour, the method offers a pragmatic complement to traditional design-oriented analyses.

4.4. Numerical Robustness of the Synthetic Sea Representation

A targeted mini-ablation study was conducted to assess the numerical robustness of the synthetic sea surface generation procedure. For the one-dimensional formulation, the prescribed significant wave height $H_s$ and peak period $T_p$ are reproduced with median relative errors typically below $1\%$, demonstrating strong statistical convergence of the random-phase spectral synthesis. Sensitivity analysis indicated that simulation duration primarily affects the accuracy of $H_s$, consistent with the stochastic convergence properties of variance-based statistics, while the peak enhancement factor and frequency resolution influence the stability of $T_p$ estimation.

For the two-dimensional directional extension, median relative errors in $H_s$ evaluated at a single spatial grid point are on the order of 7–$13\%$ for the baseline configuration. This deviation is physically expected due to spatial interference effects inherent in directional superposition and does not indicate spectral inconsistency. Frequency resolution is identified as the dominant numerical control parameter, whereas variations in directional discretisation and the spreading parameter produced comparatively minor changes in bulk statistics.

These findings confirm that the selected synthesis parameters provide a stable compromise between computational efficiency and statistical fidelity. Importantly, the ablation study demonstrates that the synthetic sea generation does not introduce systematic bias into the sea-state descriptors used to interpret motion envelopes, thereby reinforcing the methodological consistency of the proposed framework.

4.5. Limitations and Scope for Future Work

Several limitations of the present study should be acknowledged. First, the analysis relies on sea-state descriptors derived from operational monitoring data rather than from fully resolved directional wave spectra. Although wave direction is available in the dataset and used for directional analysis, harbour environments are characterised by complex wave transformation processes, including reflection, diffraction, and resonance. As a result, the effective wave-forcing direction acting on a berthed vessel may differ from the incident offshore wave direction. Consequently, directional trends identified in this study should be interpreted as statistical relationships rather than deterministic representations of the local wave field.

Second, each motion degree of freedom is analysed independently. This approach maximises the use of heterogeneous field data but does not explicitly account for coupling between motions such as sway–yaw or roll–pitch interactions. In practice, moored vessel responses are governed by coupled hydrodynamic and mooring dynamics, and these interactions may influence the joint behaviour of the motions under certain sea-state conditions.

Third, the response envelopes represent aggregated behaviour across multiple vessels and berth locations. Vessel-specific parameters such as hull geometry, loading condition, and mooring configuration are not explicitly normalised in the analysis. As a result, the derived envelopes should be interpreted as representative statistical characterisations of the response rather than vessel-specific predictions.

Future work could extend the proposed framework in several directions. Incorporating directional wave spectra or higher-resolution wave measurements would improve characterisation of the directional forcing and enable separation of wind-sea and swell components. Further developments could also include vessel classification strategies or normalisation approaches that account for differences in ship size, loading condition, and mooring configuration. Finally, probabilistic or multivariate modelling approaches could be explored to capture coupling between motion components while preserving the data-driven philosophy of the present framework.

5. Conclusions

This study developed a statistical framework for analysing ship motions at berth using heterogeneous long-term field measurements collected across multiple vessels and berth locations. The proposed approach enables sea-state-conditioned characterisation of berth motions without requiring simultaneous measurements across all motion degrees of freedom or detailed vessel-specific hydrodynamic models.

A robust two-stage quality-control procedure combining physical plausibility screening and regression-based outlier detection proved effective in removing inconsistent response–sea-state pairs while preserving physically meaningful high-response observations. This filtering step significantly improved the coherence of the response–sea-state relationships and provided a reliable basis for subsequent statistical analysis.

Sea-state-conditioned response envelopes derived from the quality-controlled dataset revealed clear relationships between wave severity and berth motion behaviour. Motion variability increased systematically with significant wave height across all degrees of freedom, confirming wave energy as the primary driver of berth motions. Quantitative envelope metrics further showed that the surge exhibits the strongest translational sensitivity to wave height. In contrast, roll displays the largest normalised response coefficient ($C_m\approx 0.39$), indicating strong amplification relative to wave height. Rotational motions, particularly roll and yaw, also exhibited the largest variability between typical and extreme responses, highlighting their susceptibility to intermittent amplification.

Quadratic envelope analysis demonstrated that several motion responses exhibit nonlinear growth as wave height increases. In particular, roll, yaw, and surge showed clear superlinear behaviour in the median response envelopes, whereas extreme responses for heave and sway exhibited stronger curvature in the upper-percentile envelopes. These findings emphasise the importance of distinguishing between typical and extreme response behaviour when assessing berth operability.

Directional analysis of measured wave directions reveals systematic directional dependence across several motion components, particularly for horizontal and rotational motions. This behaviour reflects the influence of berth orientation, mooring configuration, and local wave transformation processes within harbour environments.

Synthetic sea surface elevations generated using a spectral random-phase approach provided statistically consistent time-domain representations of the analysed sea states. Validation confirmed that the generated surfaces reproduced the prescribed spectral properties and sea-state parameters, enabling physically interpretable wave realisations linked to the derived response envelopes.

Overall, the proposed methodology demonstrates that statistically robust, field-based analysis can provide valuable insight into berth motion behaviour even when datasets are incomplete or heterogeneous. By avoiding restrictive assumptions about motion simultaneity and vessel-specific modelling, the framework provides a transparent and practical basis for screening-level berth operability assessments and comparative analyses across ports and vessels.

Future work may extend the framework by incorporating directional wave spectra, vessel classification or normalisation strategies, and probabilistic models to capture coupling between motion components. Such developments would further strengthen the integration of data-driven berth monitoring with traditional hydrodynamic analysis approaches.

Appendix A.1. Statistical Interpretation of Seed Variability

The synthetic sea surface is generated as a finite random-phase superposition of harmonic components. Although the target spectrum is normalised such that the theoretical variance satisfies $m_0=H_s^2/16$, the variance estimated from a finite-duration realisation remains a random variable.

For a Gaussian process, it is well known from classical statistical theory that the scaled sample variance follows a chi-squared distribution. Accounting for correlation in the time series through an effective number of independent degrees of freedom ${\nu }_{\mathrm{eff}}$, the variance estimator approximately satisfies

$${\displaystyle \frac{{\nu }_{\mathrm{eff}}{\widehat{\sigma }}_{\eta }^2}{{\sigma }_{\eta }^2}}\sim {\chi }_{{\nu }_{\mathrm{eff}}}^2,$$

(A1)

which implies

$$\mathrm{Var}\left({\widehat{\sigma }}_{\eta }^2\right)\approx {\displaystyle \frac{2{\sigma }_{\eta }^4}{{\nu }_{\mathrm{eff}}}}.$$

(A2)

Since $H_s=4\sqrt{{\sigma }_{\eta }^2}$, first-order error propagation yields

$${ϵ}_{H_s}\approx {\displaystyle \frac{1}{2}}{\displaystyle \frac{{\widehat{\sigma }}_{\eta }^2-{\sigma }_{\eta }^2}{{\sigma }_{\eta }^2}}=\mathcal{O}\left({\nu }_{\mathrm{eff}}^{-1/2}\right).$$

(A3)

For the 2D synthesis, variance is evaluated at a single spatial grid point over a relatively short duration (600 s), leading to a smaller ${\nu }_{\mathrm{eff}}$ compared to the 1D one-hour synthesis. Consequently, fluctuations of approximately 7–13% in reconstructed $H_s$ across independent phase realisations are statistically expected and reflect finite-sample Monte Carlo variability rather than systematic bias in spectral normalisation.