A Bayesian inverse dynamic approach for impulsive wave loading reconstruction: Theory, laboratory and field application

The measurement of wave forces acting on marine structures is a complicated task, both during physical ex- periments and, even more so, in the field. Force transducers adopted in laboratory experiments require a minimum level of structural movement, thus violating the main assumption of fully rigid structure and introducing a dynamic response of the system. Sometimes the induced vibrations are so intense that they completely nullify the reliability of the experiments. On-site, it is even more complex, since there are no force transducers of the size and capacity able to measure such massive force intensity acting over the very large domain of a marine structure. To this end, this investigation proposes a Bayesian methodology aimed to remove the undesired effects from the directly (laboratory applications) or indirectly (field applications) measured wave forces. The paper presents three applications of the method: i) a theoretical application on a synthetic signal for which MATLAB ® procedures are provided, ii) an experimental application on laboratory data collected during experiments aimed to model broken wave loading on a cylinder upon a shoal and iii) a field application designed to reconstruct the wave force that generated recorded vibrations on the Wolf Rock lighthouse during Hurricane Ophelia. The proposed methodology allows the inclusion of existing information on breaking and broken wave forces through the process-based informative prior distributions, while it also provides the formal framework for uncertainty quantification of the results through the posterior distribution. Notable findings are that the broken wave loading shows similar features for both laboratory and field data. The load time series is characterised by an initial impulsive component constituted by two peaks and followed by a delayed smoother one. The first two peaks are due to the initial impact of the aerated front and to the sudden deceleration of the falling water mass previously upward accelerated by the initial impact. The third, less intense peak, is due to the interaction between the cylinder and remaining water mass carried by the individual wave. Finally, the method allows to properly identify the length of the impulsive loading component. The impli- cations of this length on the use of the impulse theory for the assessment or design of marine structures are discussed.

The measurement of wave forces acting on marine structures is a complicated task, both during physical experiments and, even more so, in the field. Force transducers adopted in laboratory experiments require a minimum level of structural movement, thus violating the main assumption of fully rigid structure and introducing a dynamic response of the system. Sometimes the induced vibrations are so intense that they completely nullify the reliability of the experiments. On-site, it is even more complex, since there are no force transducers of the size and capacity able to measure such massive force intensity acting over the very large domain of a marine structure. To this end, this investigation proposes a Bayesian methodology aimed to remove the undesired effects from the directly (laboratory applications) or indirectly (field applications) measured wave forces. The paper presents three applications of the method: i) a theoretical application on a synthetic signal for which MATLAB® procedures are provided, ii) an experimental application on laboratory data collected during experiments aimed to model broken wave loading on a cylinder upon a shoal and iii) a field application designed to reconstruct the wave force that generated recorded vibrations on the Wolf Rock lighthouse during Hurricane Ophelia. The proposed methodology allows the inclusion of existing information on breaking and broken wave forces through the process-based informative prior distributions, while it also provides the formal framework for uncertainty quantification of the results through the posterior distribution.
Notable findings are that the broken wave loading shows similar features for both laboratory and field data. The load time series is characterised by an initial impulsive component constituted by two peaks and followed by a delayed smoother one. The first two peaks are due to the initial impact of the aerated front and to the sudden deceleration of the falling water mass previously upward accelerated by the initial impact. The third, less intense peak, is due to the interaction between the cylinder and remaining water mass carried by the individual wave.
Finally, the method allows to properly identify the length of the impulsive loading component. The implications of this length on the use of the impulse theory for the assessment or design of marine structures are discussed.

Introduction
Impulsive loading due to a breaking wave or to the initial impact of a broken wave is of great interest for the design of offshore and coastal structures. The transient nature of this load, relatively short duration (e. g. some 0.02 s (Goda et al., 1966)) and high intensity, makes it of great interest not only from the hydraulic point of view but also from the structural one, Dermentzoglou et al. (2020). The time domain representation of impulsive loading is characterised by sharp shapes that are not adequate to properly highlight its particular nature and dangerousness. However, a frequency domain approach better serves to present how the content of energy within an impulsive load can be dangerous for every kind of structure. Indeed, the energy is spread among a large range of frequencies (theoretically from 0 to ∞, e.g. Fig. 1) so that the risk for induced resonance, and consequently, amplification of the effective load, is significant. Fig. 1 shows, through the well-known Kronecker delta, these phenomena wherein the upper panel presents the time series of the Kronecker delta and the lower panel its Hilbert-Huang spectrum (HHs) with energy spreading through all the analysed frequencies.
This is true not only for on-site conditions, where the amplification of the effective load can be detrimental for the structure integrity (Serinaldi and Cuomo, 2011), but also when the impulsive wave force has to be measured in a hydraulic laboratory. As a result, the measurement of wave forces acting on marine structures is a complicated task, both during physical experiments and, even more so, in the field. Force transducers adopted in laboratory experiments require a minimum level of structural movement, thus violating the main assumption of fully rigid structure and introducing a dynamic response of the system that masks the hydrodynamic load (Dassanayake et al., 2019a). On-site, it is even more complex, since there are no force transducers of the size and capacity able to measure such massive force intensity acting over the very large domain of a marine structure. Field wave pressures have been measured to determine the overall loading, producing benchmark information for understanding the interaction between wave and structures (Bullock et al., 2007). Wave pressures have been measured with success in several experimental campaigns (Cuomo et al., 2010;Cuomo et al., 2007;de Almeida and Hofland, 2020;de Almeida et al., 2019;Stagonas et al., 2016); however, the overall description of the total wave forces is affected by several assumptions and possible inaccuracies in the List of symbols d (t) displacement of the mass or measured force h (t− τ) or IRF unit-impulse response functions F (τ) external system perturbation M mass (or equivalent mass) of the modelled body k dimensional system stiffness c the dimensional system viscous damping coefficient ω n system natural frequency Toeplitz matrix representing the convolution operation when the exerted impact hammer force is known and the system IRF is unknown R matrix resulting from the QR decomposition H m0 spectral significant wave height H m0MAX maximum recorded significant wave height T P spectral peak wave period T S spectral significant wave period, T S = T P /1.07 T m = spectral mean wave period T S = T P /1.19 d P spectral peak wave direction CH i accelerometer ith signal H 0. 21% wave height with exceedance probability equal to 0.21% spatio-temporal integration of those pressures Martinelli et al., 2018). Therefore a force reconstruction method is required with its inherent downsides due to the solution of the underlying inverse problem (Maes et al., 2018;Sanchez and Benaroya, 2014). The problem of characterising the impulsive wave loadings has attracted researchers' interest since 1958, when Hall performed the first laboratory experiments aimed to characterize the breaking wave forces on a circular pile located on a sloping beach, (Hall, 1958). Several authors follow these pioneering tests and are nicely summarized by Tu (2018) and Tu et al. (2017a). Goda et al. (1966) were the first to formalize a mathematical model to describe breaking wave loading and to highlight the need to consider within the description of the loading condition the dynamic of the affected structure, at both the model and prototype scale. Goda based his final formula on the combination of experimental results and von Karman theory (Von Karman, 1929). Later (Campbell, 1980), performed drop tests instead of wave impact tests, in order to achieve a reasonably large Froude number (e.g. > 0.6) so that the total load was mainly dominated by the slamming component. Despite achieving high rigidity in test set-up (natural frequency around 550 Hz), the amplification due to the dynamic response was large enough to mask the hydrodynamic load. Hence, to properly describe the pure hydrodynamic load, the experimental system was modelled as a lumped mass and was forced with a hyperbolic function. The goal of the method was to identify which shape of the hyperbolic function allowed the best match between the dynamic response of the experimental set-up and the response of the single degree of freedom (SDoF) model. More recently (Wienke and Oumeraci, 2005), performed a large scale test aimed at identifying the breaking wave loading on vertical cylinder under the action of focused wave groups. The model comprised a cylinder installed in deep water and fastened at both ends. Also in this case the transient nature of the impulsive wave loading induced dynamic response of the experimental set-up, so they applied a similar method to that of Campbell (1980), though removing a known quasi-static force from the experimental dynamic response. The approach allowed the verification of the assumed theoretical description of the impulsive load by contrasting a SDoF model response and the measured dynamic force. However, both Campbell (1980) and Wienke and Oumeraci (2005) avoided the use of the inverse method, and under the hypothesis of linear response of the experimental structures, instead verified their formula by means of the convolution process between the developed empirical equation and the impulse response function (IRF) of the structure. Dynamic amplification of breaking wave loading during experimental study was also highlighted by Choi et al. (2015), where they quantified, and subsequently removed, this unwanted effect by means of the combined use of the Empirical Mode Decomposition (EMD) and Computational Fluid Dynamic model (CFD). Later, a similar approach based on Ensemble Empirical Mode Decomposition (EEMD) was applied by Dassanayake et al. (2019a) and Dassanayake et al. (2019b) to remove the effects of the vibration induced by broken waves on the experimental set-up aimed to model an offshore rock lighthouse. Despite the EEMD approach being more accurate in removing the dynamic response of the structure than EMD, it still presents disadvantages related to the overestimation of the quasi-static force component, as also highlighted by Tu (2018).
It was only within the WaveSlam project (HYDRALAB IV framework) that the inverse estimation of the breaking wave force acting on marine structures has been successfully undertaken. Four different approaches were proposed within Tu's related PhD thesis (Tu, 2018), all of them based on the deconvolution between the recorded wave force and the dynamic response of the structure. The first method, called optimization-based deconvolution (OBD) (Tu et al., 2015), is based on the minimization of the Euclidean distance between the measured force and the modelled force given by the superimposition of the scaled and shifted hammer test responses. The other three methods are mainly based on the deconvolution between the IRF and the measured dynamic response of the structure under the breaking wave action. The so-called "horizontal approach" is a two-step approach based on the impact hammer test results. The IRF is identified by means of conjugate gradient technique while it is later applied to the recorded dynamic response through a weighted eigenvector expansion method in order to obtain the wave slamming force. The horizontal approach relies on two main parameters for the regularization of the solution, the stopping factor and the weighting factor for the first and second steps respectively. Both are defined by the user in order to control the noise effect in the IRF definition (first step) and to discriminate the smallest eigenvalues of the deconvolution matrix (second step), thus reducing the risk of numerical instability through a regularization approach (Tu et al., 2017b). The so-called "vertical approach" uses the linear regression technique. Similarly to the horizontal approach, it also reconstructs the wave impact force at each investigated location by using the hammer impact force, the hammer response force, and the wave response force at the same measurement location. For each investigated location, the wave impact force is conceived of as a result of the hammer hitting this location with different amplitudes many times in a row, hence the interval between every two imaginary hammer impacts is an input required from the user and is called step factor. The step factor indirectly controls the size of the deconvolution matrix and the accuracy of the reconstructed impact force. The last proposed approach is the "extended vertical approach". Similarly to the OBD, the extended vertical approach accounts for the contribution of the impacts at different locations into the measured force responses, while treating each transducer simultaneously. In this approach, the response locations and the impact locations are distinguished. More recently Maes et al. (2018) applied a recursive joint input-state estimation algorithm for the inverse estimation of the breaking wave loading on hydro-elastic model scaled wind turbine monopile and the induced members forces at the base of the flexible structure. The algorithm is based on the dynamic behaviour of the flexible monopile along the incoming wave direction where the modal parameters are experimentally identified via impact hammer tests. The results show relatively close agreement between the measured and reconstructed forces with an average absolute error around 27% for the impact force and 19% for the overturning moment. However, the overall method relies on the assumption of triangular pressure distribution for which there is no evidence that it can be used within the whole loading process. Despite this assumption, Maes et al.' work sets a foundation for the inverse wave force identification through the dynamic response of the structures. Finally, to estimate the magnitude of the slamming load on offshore wind turbine, Paulsen et al. (2019) applied a simplified dynamic model of the laboratory set-up to describe the transfer function and hence partially remove the dynamic oscillation of the cylinder. However, the methodology is not described in detail by Paulsen et al. because of the different focus of the analysis.
This work intends to make progress in the application of the inverse method to reconstruct wave forces exerted on marine structures. We aim to present a new Bayesian inverse method to reconstruct both field and laboratory forces due to breaking or broken waves. While tackling the three main downsides of the inverse methods, i.e. solution existence, uniqueness and stability (Aster et al., 2018), the proposed approach will provide not only a proper framework to analyse future laboratory and field data from offshore and coastal structures, but also a tool to account for the prior knowledge on breaking wave forces and a formal approach for the uncertainty quantification of the results through the posterior distribution. Therefore, this paper is not aimed at producing a comprehensive description of the specific impulsive load due to the broken waves on a cylinder upon a shoal, but instead at presenting and describing a useful Bayesian methodology to achieve a more comprehensive and general result. In order to achieve this, the paper presents an introductory overview on the experimental problems related to the measurement of the wave forces and on the issues connected with the inverse methodology required to solve the inevitable violation of fully rigid model assumption, based on linear systems theory. The convolution between the input signal and the IRF is the core concept of the methodology, therefore, the assumption of linear elastic behaviour of the structure is implicitly adopted. From the statistical point of view, the main hypothesis is related to the noise affecting the data that is considered normally distributed and independent. Through the development of the paper, the method is applied to both laboratory and field data. Despite the different nature of the recorded experimental and field vibrations -one is a laboratory effect whereas the other is a real structural response -the methodology is successfully applied to reconstruct both wave forces.
The method will be presented in the following chapter where the main theoretical background and numerical issues will be treated in order to provide the required background knowledge. The proposed method will then be applied to the laboratory force measurements (chapter 4) and field accelerations measurement (chapter 5) as illustrative examples of application. Finally, discussions about the main results of the applications and the resolved issues will be gathered in chapter 6.

Method
The proposed solution for the inverse problem is a merger of structural and statistical models, thus it requires a proper formulation of the fundamental hypotheses for both aspects. The approach describes the dynamic behaviour of the investigated structure by means of damped single degree of freedom model (SDoF) (whether it is laboratory or prototype one), under the main structural hypothesis of linear timeinvariant behaviour. This means that the wave loading that is to be reconstructed, cannot generate plastic deformation and also that the structure has fully elastic behaviour under such loading. This allows the calculation of the response (d (t) , e.g. displacement of the mass or measured force) to an arbitrary time-varying external perturbation (e.g. the wave force) by means of the superposition of a series of unit-impulse response functions (IRF or h (t− τ) ) due to a series impulses composing the external perturbation (F (τ) ). This concept is well known within the earthquake engineering as Duhamel's integral (Rajasekaran, 2009) or more generally under the mathematical concept of Fredholm integral equation of the first kind (Aster et al., 2018) and is represented in eq. (1): while the displacement IRF for damped SDoF can be written as shown in eq. (2).
where M is the mass (or equivalent mass) of the modelled body, k and c are the dimensional stiffness and viscous damping coefficient, ω n is the natural frequency calculated through the well-known equation ω n = √ . Most often in civil engineering the damping ratio does not exceed 20%, thus the damped and natural frequencies tend to be the same (Lee et al., 2018;Martinelli and Lamberti, 2011;Rajasekaran, 2009). However, despite the optimum laboratory set-up making use of stiff instruments, the nature of the connection of the sought instruments with additional elements (as in the case presented in this paper) could perhaps introduce damping to a level requiring distinction of damped and natural frequencies.
The proposed method aims to solve the inverse deconvolution operation that will remove the dynamic effect of the structure (i.e. h (t− τ) in eq. (1)) and allows the reconstruction of the wave force F (τ) . As an illustrative example, we can assume that the measured data (d (t) ) is the response of a force transducer to an external load F (τ) that we want to reconstruct by removing the dynamic response due to the model set-up. Standard laboratory force measurements rely on transducers that integrate strain gauges, thus the structure must be free to move, hence violating the hydraulic modelling assumption of a fully rigid structure. Therefore, to reconstruct a force resulting from the dynamic response of the system, e.g. a force transducer connected with a structure, the required IRF should not be expressed in term of displacement per unitary impulse (i.e. m/Ns) as in eq. (2), but in term of inertia force per unitary impulse, (i.e. N/Ns), Fig. 2. Knowing the analytical expression of the displacement IRF, the calculation of the force IRF is easily achievable by means of the multiplication between the laboratory structure mass (M) and the second time derivative of eq. (2). Eq. (3) shows the resulting expression, in which M can be simplified as presented in eq. (5). For the illustrative example, the numerical values of the dynamic parameters are described in Fig. 2.
The convolution integral in eq. (1) can be rewritten in more convenient matrix notation as presented in eq. (4): where the symbols in bold denote a vector of values varying in time, i.e. d is the recorded force response time series, m is the sought unknown external load and G is a square Toeplitz matrix representing the convolution operation. G comprises lagged IRFs, so that the rows are time-reversed and the columns are non-time-reversed versions of the IRF lagged by i and j as shown in eq. (5). The resulting convolution matrix with example columns are presented in Fig. 3.
The fundamental statistical hypothesis of the method is related to the noise affecting the data (d (t) ), that is assumed to be normally and independently distributed with a corresponding diagonal covariance matrix C D . Moreover, the proposed approach relies on Bayes' theorem as presented in eq. (6). Hence, the sought model is assumed to be a random variable so that the final solution is a probability distribution q(m|d) for the model parameters, often called the posterior distributions.
where p(m) denotes the prior distributions and f(d|m) is the conditional probability, that, given a particular model m, corresponding data d will be observed. In other words, and more specifically related to the inverse force reconstruction, we aim to identify a probability distribution for each instant (t i ) described by the time vector (t) that our identified model (m (ti) ) might have generated the measured response (d (ti) ). Furthermore, we want to probabilistically describe how our model (m (ti) ) is effective in modelling the real unknown wave force (F (τi) ) that generated the measured response d (ti) . The Bayesian approach allows the natural incorporation of the prior information about the final solution that comes from previous knowledge or experience by means of the timevarying prior distributions p(m (ti) ). Therefore, the previously devel-  ) can be directly considered within the analysis. The prior distributions are assumed to be normal distributions varying with time as in eq. (7), with expected value m prior and associated covariance matrix C M .
The likelihood that given a particular model, a response vector (d) will be observed is expressed by the likelihood function f(d|m) eq. (8).

Fig. 2.
Analytical force IRF for a linear system having mass (M) equal to 20 kg, natural (ωn) and damped (ωd) frequency the same and equal to 10 Hz, stiffness (k) equal to 1000 Nm and damping ratio (ζ) 2%, sample time step (Δt) equal to 0.001 s.
Therefore, through the resolution of the integral in eq. (6) the prior and posterior distribution are related in a way that makes the computation of q(m|d) possible. The final result of the method is a series of normal distributions (i.e. one for each sampled data value), describing the unknown wave force as shown in eq. (9).
It must be noted that q(m|d) does not provide a single value that we can consider "the wave force", thus to provide a single model output of the wave force, the maximum a posteriori value (MAP), i.e. the wave force associated with the largest value of q(m|d), is proposed as suggested by Aster et al. (2018), leading to a simplification of eq. (9) in eq. (10) and eq. (11), as proposed by Tarantola (1987).
Rewriting C − 1 M and C − 1 D in terms of matrix square root by means of, for example the Singular Value Decomposition (SVD), the MAP solution can be now calculated by the minimization of the exponent in eq. (9) resulting in a standard linear least-squares problem presented in eq. (12). min and m can be seen as a transformation that makes the data (affected by random noise) and the unknown models (intrinsically stochastic due to the Bayesian nature of the methodology) independent with a normalised standard deviation for both the data and model space respectively.

Theoretical example
Usually the convolution matrix as given by eq. (5) and shown in Fig. 3 is mildly to severely ill-conditioned, hence the inverse problem is  not straightforward because we can anticipate a severe amplification of the noise contained within the real data. In the example proposed above the condition number of the matrix G is slightly larger than 88 ′ 800 due to its wide range of singular values between 24 and 2.7 × 10 − 4 , Fig. 4. Therefore, even assuming optimistically the recording signal is affected by a 0.01% noise level, the results of the inverse process will be dominated by the amplified noise.
For explanation purposes, we can conceive the reconstruction of an impulsive wave force (black line in Fig. 5) acting on the laboratory setup characterised by the force IRF in Fig. 2 and that the force measurement (i.e. the system response) is affected by some white noise (blue line in Fig. 5). A standard approach to tackle this inverse problem would be through the application of the least-squares method with the support of the Singular Value Decomposition (SVD). Indeed, the system presented in eq. (4) (i.e. d = Gm) can be solved for m once the inverse of the matrix G is obtained via the SVD decomposition as presented in eq. (13).
where V and U are the right and left singular vector matrices respectively and S is the singular values matrix. The reconstruction of the impulsive wave force from the theoretical measurement (i.e. the noiseless red dotted line in Fig. 5) fits its noiseless data perfectly, being essentially identical to the original impulsive force (Fig. 6a). If the same procedure is applied to more realistic white noise data (i.e. the blue line in Fig. 5, C D = 0.65⋅I) the solution is meaningless. The information about the original impulsive wave force is overwhelmed by the noise, enormously amplified by the inversion process (Fig. 6b).
To control the unstable character of the proposed inverse problem, a first preliminary method can rely on the property of the Fourier transform. Indeed, the Fourier transform of a convolution between two elements is equal to the product of the two Fourier transforms, so that the solution is trivial within the frequency domain. However, the solution of the inverse problem remains extremely sensitive to small changes in the records (d) and requires a regularization process that can be achieved by imposing equal to zero the smallest elements of the Fourier transform of the records, hence obtaining a sort of truncated Fourier transform. Even  though this method is rather effective in term of computational time, does not come without downside aspects. The threshold that defines the level of the "smallness" of the Fourier transform elements to be zeroed is unknown and depends on the noise realisation, therefore for each case it should be properly defined introducing a subjective selection of this fundamental parameter. The result of this simplified method is presented in Fig. 6.c with the cyan colour.
To overcome the instability issue due to the ill-conditioning of the convolution matrix a truncated SVD (i.e. the inverse of the matrix G is obtained by using only the largest singular values) can be applied in order to reconstruct a reasonable estimation of the original force and avoid the subjective selection of the above mentioned threshold. Fig. 6.c shows the reconstructed force obtained using the L-curve criterion (Hansen, 1992;Hansen, 2007) as a guide for selecting the Tikhonov regularization parameter (Tikhonov and Goncharsky, 1987) where only the first 168 singular vectors of the matrix V are used within the inverse process. The SVD truncated approach allows the detection of the essential features of the original impulsive force; however, this technique, as well as the truncated Fourier transform, introduces some spurious oscillations and loss of resolution generating a wider impulse and reduced amplitude as shown in the zoomed box in Fig. 6.c.

Informative prior distributions
In order to properly apply the previously developed Bayesian method, the prior distributions need to be defined. Having described the dynamics of the system, this information can be incorporated into the process by means of informative priors. It is reasonable to believe that the force shown in Fig. 5 should have been applied to the system at least some instants before the change of its status (i.e. t = 0.5 s). After that, it is reasonable to assume that the maximum force value (calculated, for example, using the Wienke and Oumeraci (2005) approach) should have been applied to the structure at least a short time after the maximum response value (i.e. t = 0.57 s) after which the incident force should have dropped to 0. Moreover, assuming that the incident force is an impulsive wave force there is enough knowledge (e.g. Cuomo et al. (2010); Cuomo et al. (2007); Wienke and Oumeraci (2005)) to believe that the rising slope is steeper than the decreasing one. Hence we can assume a trapezoidal-shaped zero-covariance prior distribution that preferentially concentrates the model structure around the instant of the maximum response by imposing a zero prior with small standard deviation far from the maximum system response (Fig. 7). Fig. 8 shows the comparison of the Bayesian inverse approach and the least square method approach. It is evident that the Bayesian solution is still severely affected by some noisy oscillation and large uncertainties around the reconstructed force. Because the prior distribution has zero covariance, the resulting model realisations are quite rough.
Therefore, the prior distributions can be designed to enforce a smoothness constraint on the realisation of the posterior distribution by specifying a non-diagonal prior correlation matrix (C R ). The positive definiteness of C R can be guaranteed by constructing the matrix so that each column is the autocorrelation of a pre-selected function in which the zero-lag (unit) maximum is centred on the diagonal of C R . Due to the particular shape of the impulsive wave loading, i.e. close to a triangle with the highest corner corresponding to the peak force, the autocorrelation of a triangle function that produces a cubic approximation to a Gaussian function is selected. Moreover, we can base the correlation time scale on the previously developed knowledge of the impulsive wave force duration (Cuomo et al., 2010;Goda et al., 1966;Wienke and Oumeraci, 2005) so that the prior correlation function falls off with a time scale of 0.05 s, i.e. the zero-lag (unit) maximum of the correlation sequence is centred on the element i and zero at approximately i ± 0.05 s (Fig. 9). Hence, given the non-uniform diagonal elements of the covariance matrix C M and the correlation matrix C R , the non-diagonal prior covariance matrix C M,R is defined as in eq. (14):    10 shows the final result of the developed Bayesian inverse method, in which full use of the previous knowledge about both the dynamic behaviour of the system and the physical knowledge about the breaking wave loading have been incorporated within the prior distributions and therefore within the inverse process. The obtained solution is, not surprisingly, considerably improved by the more restrictive prior model. This is because the true model is highly restricted and thus consistent with the prior distributions. Moreover, it can be recognised that, despite the restrictive prior model, the information carried by the data is not overwhelmed by the prior distribution. The zoomed box in Fig. 10 clearly depicts a smaller credible interval for the rising part of the impulsive force (i.e. 0.5 ≤ t ≤ 0.538) than for the falling one (i.e. 0.38 < t ≤ 0.8) even if the assumed prior variance is the same. Thus, it can be argued that the slope of the original signal can be interpreted as an index of the relative importance or strength, within the Bayesian process, between the information carried by the data and that carried by the prior model. The 95% credible interval is not the standard 95% confidence interval, rather it is the 95% probability interval calculated from the posterior distributions, so that there is 95% probability that each m (t i ) value lies within the corresponding symmetric interval around the MAP value.
In this chapter, three different methods to solve the inverse problem aimed to reconstruct the incident wave force from a noisy signal recorded on a, or from a structure have been presented. The first and more simplified method makes use of the Fourier transform properties and the subjective selection of the threshold to solve the instability issue due to the noise affecting the records. The second is based on the SVD decomposition of the convolution matrix, the Tikhonov regularization criterion and the L-curve technique for the selection of the threshold aimed to identify the largest singular values to consider within the  A. Antonini et al. resolution of the deconvolution process. The third, and most comprehensive one, is based on the Bayes theorem, it provides the framework to inherently quantify the uncertainty in the final solution and to account for the prior knowledge about impulsive wave loading and structural dynamic. All the proposed methods allow a reasonable reconstruction of the incident force, each of them with its own downsides, but still all of them applicable to further case studies involving coastal and offshore structures under impulsive wave loadings.
The above theoretical example, as well as the following laboratory and field applications, have been entirely developed within the MAT-LAB® environment. The entire procedure chain adopted to complete the theoretical example is also released as additional material to this document, aiming for a straightforward easy application to different case studies.

Laboratory setup
A series of physical model tests on a vertical cylinder upon a variety of 2D shoals was performed within the STORMLAMP (STructural behaviour Of Rock Mounted Lighthouses At the Mercy of imPulsive waves) research project framework in the wave flume of the COAST Laboratory, University of Plymouth. The flume is 35 m long, 0.6 m wide and 1.2 m high. A lighthouse is modelled as a vertical aluminium cylinder (weight 9.88 kg) installed at the middle of the shoal, while the adopted foreshore steepness for the test reported here is 1:5. The water surface is measured by means of 16 wave gauges (WG), spread before and after the shoal, while two cameras, standard and high-speed ones, are used to capture the wave development along the foreshore (standard ones) and at the cylinder (high speed one) ( Fig. 11 and Fig. 12). The high speed camera records were also used to evaluate the runup along the cylinder by means of the methodology presented by Dassanayake et al. (2019a). The offshore flume bed is flat and the mean free surface is coincident with the upper part of the shoal, i.e. 0.5 m (Fig. 11). The cylinder diameter (D) is 0.12 m, while the width of the upper shoal platform is 0.36 m. The 0.5 m high cylinder is suspended from its top and behaves as a vertical cantilever, leaving a minimal gap (i.e. 0.7 mm) between the cylinder bottom surface and the shoal. The top of the cylinder is connected to a 6 degrees of freedom load cell (model:6A40B-500/20 -weight 0.4 kg (Interface, 2019),), that in turn is connected to a beam which is part of the main supporting structure (Fig. 12). The height of the force transducer is 40 mm, while the origin of the coordinate system is located 32 mm above the cylinder top surface, so that the cylinder bottom surface is at 532 mm from the origin. The load cell is equipped with 6 temperature-compensated bridges providing output for each of the 6 degree of freedom. Therefore the output signals have to be post-processed by means of a 6 × 6 calibration matrix in order to extract the force and moment values. The set-up enables force measurements along three perpendicular axes with three simultaneous moments. The sample frequency is 5120 Hz, but all the signals have been decimated to 1000 Hz in order to reduce the computational effort of the inverse process. Regular, irregular and focused waves were run; however, this investigation considers the applied methodology to reconstruct the wave force, hence results from regular waves only are presented. The assumption of 2D model is valid for the present test, hence, the wave force is acting along the negative y-direction and the induced moment is positive around the x-direction, Fig. 11.
The cantilever scheme leads to a versatile set-up but also to an unavoidable reduction in the overall system stiffness, requiring the need to properly address the wave-induced vibrations. Fig. 13 upper and middle panels show the typical recorded force and moment patterns for a regular wave case characterised by wave height (H) equal to 0.14 m measured at WG11 on the shoal (water depth 0.1 m and distance from the cylinder 0.5 m) and period (T) equal to 1.5 s (Fig. 13 dotted line  lower panel). A similar vibration pattern is present for all the wave states, highlighting the overwhelming effects of the structural dynamic response on the recorded force. Fig. 13 also shows the raw records and the Hilbert-Huang spectrum (HHs) presenting a clear pattern due to the natural frequency of the laboratory set-up which becomes the dominant feature of the records. Between 3.0 and 3.2 s a sudden jump in the instantaneous frequency and energy is detectable which is likely to indicate the instant at which the wave impacted the structure, as confirmed by the measured runup (lower panel solid line). Proceeding along the signal development, the natural frequency of the structure becomes dominant as shown by the instantaneous energy concentration between 12 and 12.5 Hz. Less energetic intrinsic mode functions (IMFs) are also grouped around the lower frequencies close to the incoming wave frequency equal to 0.66 Hz, however due to the large difference in instantaneous energy they are barely discernible. Overall, it is evident that the recorded force is not the wave force but the response of the model to an external perturbation. Finally, we want to stress that, although in this example the dynamic response is particularly effective in corrupting the measurement, any impulsive wave force measurements should be properly post-processed with different deconvolution techniques in order to extract real features and intensity of the incident force, regardless of whether the laboratory model is relatively stiff.

Laboratory Bayesian inverse method application
As for any dynamic system, the first step is the identification of the dynamic behaviour. In order to properly describe the dynamic response of the laboratory model, impact hammer tests have been performed with the aim to experimentally reconstruct the force IRF. The impact hammer tests made use of a piezoelectric impact hammer equipped with a rubber head (Fig. 14) which was used to hit the dry cylinder, i.e. without any surrounding water, 3 times around a lower location where the wave impact is expected (Fig. 15). By using the dry IRF within the inverse process, we are implicitly assuming that the dynamic parameters of the laboratory model remain the same during the interaction with the wave. From the preliminary results of wet IRFs, we identified that the damping ratio increases and the natural frequency decreases due to the additional viscous damping and added mass due to the surrounding water. However, the uncertainty in the level of the water that should have been considered to properly reproduce the wave impact conditions do not allow the use of the wet IRF within the inverse process and therefore the dry IRF has been used through the entire paper. The final adopted IRF is the time average of 3 IRFs each of them calculated, as will be described, by dividing the signals shown in Fig. 15 in shorter and equally spaced segments with a length equal to 2.06 s and highlighted by the red dotted lines.
The sought IRF can be seen as the time domain image of the frequency response function (FRF) of the system, so that the IRF can be calculated as the inverse Fourier transform of the ratio between the Fourier transform of the system output (lower panel in Fig. 15) and the Fourier transform of the system input (upper panel in Fig. 15) or simply the time domain deconvolution of the two. Under the hypothesis of fully rigid cylinder and supporting structure the overall laboratory set-up can be approximated by a single degree of freedom system with a strict relation between the rotation at the force transducer and the displacement at the tip of the cylinder so that the same IRF can be used for both force and moment. Appendix A presents the derivation of the IRF and the comparison between the IRF, calculated under the above-mentioned hypothesis and by explicitly taking into account both the force and the moment. However, despite the robust theoretical basis, the operation as described above in the case of noisy discrete measurements is illconditioned, so that a regularization procedure needs to be applied also at this stage. Here, the issue is tackled by means of the least-squares solution supported by the QR decomposition. Each of the three experimental IRFs, calculated by means of the three signal segments identified in Fig. 15, is calculated as the solution of the linear system presented in eq. (15): Fig. 14. Impact hammer tests and impact location.
where h r is one of the three segments representing the response of the laboratory model to the impulsive force exerted by the impact hammer (i.e. h h ) that in turn is divided by the time integral of the impulsive force recorded by the impact hammer, i.e. h h . G h is a matrix defined using the same method of matrix G (eq. (4)) with the main difference that, in this case, each column is defined as a lagged hammer force, i.e. h h . Therefore, the sought solution is the IRF that minimises both the norms in eq. .
The applied regularization approach is aimed at treating the smallest elements along the diagonal of the matrix R (obtained from the QR decomposition of the matrix G h ) as zeros, so that the effect of the noise in the impact hammer records does not play a major role in the final solution. Plotting the elements along R's diagonal is enough to identify a reasonable regularization threshold that in this example is set to 0.5; however this value should be evaluated for each case. Fig. 16 shows the calculated average IRF that has been adopted for all the following analysis. The clear presence of multiple components with their own frequency (13, 90, 475 and 535 Hz, see Appendix A) might be related to the quasi-rigid rotation of the cylinder around the transducer (13 Hz), to the second flexural mode of the cylinder, having the centre of mass moving in phase opposition relative to the tip (90 Hz), to the first natural mode of a cylindrical cantilever element (475 Hz) and to the supporting structure vibrations (535 Hz).  Finally, in order to identify the intrinsic noise within the data, and then define the covariance matrix (C D ) associated with the data, the signal is assumed to comprise a smoothly varying function plus additive Gaussian noise with zero mean and variance to be estimated; the methodology described by D'Errico (2007) is applied here, for the estimation of the signal noise variance. The identified variance values are 0.0015 N 2 and 2.5 × 10 − 6 (Nm) 2 for the force and moment, respectively, and are assigned to the elements along C D 's principal diagonals.

Laboratory informative prior distributions
As presented in the theoretical example, the definition of the prior distributions is based on the previous knowledge on the impulsive wave loading on cylindrical structures (Goda et al., 1966;Tanimoto et al., 1987;Von Karman, 1929;Wienke and Oumeraci, 2005). Despite alternative approaches being available, recently the work of Wienke and Oumeraci has been successfully applied in preliminarily investigations of wave loading on offshore rock lighthouses (Trinh et al., 2016), hence it is used as reference for the definition of the prior distribution. However, a large proportion of the waves that interact with the lighthouse, and accordingly, also in the present laboratory experiments, rarely  break directly onto the structure, instead they mostly reach the structure already broken with an initial aerated and turbulent front (Bressan et al., 2018). Therefore a modification is applied to the standard Wienke and Oumeraci approach and the wave celerity (C b ) is calculated according to the method of Bonneton (2004) for broken waves in the surf zone, resulting in a value 1.5 m/s. Hence, the maximum value for the prior distribution is kept equal to the maximum force calculated according to the modified Wienke's method, i.e. 31 N and is approximately applied at 0.06 m from the bottom of the cylinder (i.e. 0.47 m from the origin of the axis). These dimensions are associated with a curling factor of 0.46 and wave crest of 0.08 m at the breaking point ( Fig. 13 upper panel).
Finally, due to the large uncertainties on the multitude of phenom-    ena affecting the interaction between the flow and the structure, a relatively non-restrictive standard deviation is assumed. Thus, a value equal to half the mean value is assigned to the standard deviation in order to fully describe the prior normal distributions for both horizontal force and overturning moment, as shown in Fig. 17. Due to the non-fully breaking nature of the waves action, for the definition of the prior distribution, the time length of the impulsive force is initially estimated to be 0.04 s, according to Goda et al. (1966) who presented the longest values among the available impulsive model lengths, i.e. D 2⋅C b . Hence, the prior correlation function falls off with a time scale of 0.04 s.

Laboratory results
Results of the analysis are shown in Fig. 18 for a record of 5 incident regular waves, with 1.5 s period and wave height around the breaking point of 0.14 m (a movie of the 5 incident waves and the obtained results is available as additional material to this document). Panel a) shows the water surface elevation recorded at WG11 situated on the shoal (i.e. water depth 0.1 m and distance from the cylinder 0.5 m) and the measured runup, while panels b & d) show the identified wave force and induced moment respectively, and panels c & e) their HHs. The runup is defined as the level B runup proposed by Grue and Osyka (2021), i.e. the runup of a thin layer of water and air mixture, and water layer which was no longer attached to the surface of the pile, or high spray concentration.
The dynamic amplification due to the structure is completely removed from the records allowing the description of a clear signal and the identification of the main loading features due to the broken waves. The detected force and moment highlight the presence of common features within the signals. Two peaks are clearly visible for all the loading events while a third one is slightly less pronounced but still present for all the events. The first peak is related to the violent impact of the first broken aerated front (red arrow Fig. 18b,d), the second is mainly due to the sudden deceleration of the falling down water mass previously upward accelerated by the impact with the cylinder (green arrow Fig. 18b,d). A third, less intense, peak due to the remaining water mass carried by the wave is also detected and is clearly visible for the third wave (purple arrow Fig. 18b,d) but it is also present within the other loading events with a smaller intensity. Despite the abovedescribed loading mechanism is confirmed by the movie provided with this document some doubts arise regarding the effects of the cylinder compliance and movement. The moment transducer actually measures a small transducer deformation and rotation of the cylinder that due to the moderately high frequency of the induced oscillations and distance from the hinge point may result in a significant velocity of the cylinder in contact with water. However, to what extent this process affects the reconstructed force is not trivial to define and it has been assumed negligible in this work in light of the reasonable agreement with the later-described field results. As identified by Liu et al. (2019) and Kristiansen and Faltinsen (2017) for breaking wave on a vertical deep water cylinder, the content of energy for the impulsive load part is spread over a frequency range broader than the incident wave frequency. For the analysed case the waves break before the structure, hence the front that first impacts on the cylinder is extremely turbulent inducing a longer rise time but also energy content at higher frequencies that reach up to 40 times the wave frequency (i.e. 0.66 Hz), as shown by the HHs in Fig. 18c,e. This is particularly relevant for stiff structures, like the granite masonry offshore rock lighthouses, for which the observed natural frequencies for the first two modes are in the range between 4 and 8 Hz, Brownjohn et al., 2018). Therefore, a detailed description of the higher wave load harmonics is essential to describe the induced dynamic response. Moreover, a constant low frequency component with a value close to the incident wave frequency (i. e. 0.66 Hz) is visible all along the time series shows in Fig. 18c,e. For all the impact events, a sudden jump in the instantaneous frequency is detected with a concentration of energy during the rising part of the impulsive load, particularly pronounced for the fourth event. The detected rise times range between 7% and 9% of the wave period (i.e. 1.5 s) and between 20% and 30% of the whole impulsive loading duration that lasts around 0.2 and 0.3 s. It is important to highlight this aspect because several approaches, e.g. Goda et al. (1966), (Goda et al., 1966); Wienke and Oumeraci (2005), do not consider the rise time when describing the breaking wave force time series, whereas it is indeed the part of the impulsive load where the energy is largely concentrated for the broken wave action. Fig. 19 shows details of the results for the third event in Fig. 18, in which the detected wave force application point is shown together with the runup in panel c. The runup measurement is obtained by the automated image processing method described in Dassanayake et al. (2019a), while the application point is the crude ratio between the moment and the force.
As expected, both the wave force and the induced moment are related to the runup as already highlighted by Peregrine (2003) for a vertical wall under breaking waves. The initial increase in runup (Fig. 19.c: 3.0-3.05 s) is largely due to the jet and the aerated water mass generated by the wave breaking before the cylinder, therefore little or no pressure is exerted on the structure. Subsequently, the primary front of the broken wave reaches the structure. It is projected upward by the pressure gradient due to the high pressure developed during the contact between the water mass and the cylinder (Fig. 19: 3.05-3.15 s) until it reaches the maximum runup level (Fig. 19: 3.15 s). The force application point correlates reasonably with the magnitude of the force, reaching its maximum slightly later than the maximum force, thus also inducing a different phase between the maximum force and maximum moment, then it suddenly drops as the force decreases. At the point of maximum runup, the water is in a nearly in free fall (Fig. 19: 3.15-3.25 s), exerting little pressure on the water below and resulting in the reduction of the force and moment. As it falls down, the water must be decelerated by a pressure gradient that is again supported by high pressure at the base of the cylinder and therefore by a second peak in the horizontal force and moment, (Fig. 19: 3.25-3.30 s). As expected the application point for the second peak is quite low and more steady than the first one, highlighting the non-impulsive nature of this loading cycle. The sudden drop of the runup around 3.3 s is mainly due to the inaccuracy of the video camera technique that failed to distinguish between the thin layer of water that was no longer continuously attached to the cylinder surface as shown in the additional video available with the paper. After the end of the impulsive loading component (Fig. 19: 3.30 s) the remaining part of the water mass carried by the wave reaches the cylinder generating a secondary load cycle (Fig. 19: 3.35-3.45 s) that is less violent than the primary one, but still shows a slightly impulsive nature. On average it was observed that the intensity of this secondary load cycle ranges between 15% and 25% of the primary one and it lasts for a duration that ranges between 5% and 12% of the wave period and between 40% and 60% of the primary impulsive load. However, despite the reduced intensity of the secondary load cycle, it consistently shows a relatively high application point that ranges between 60% and 80% of the runup levels; thus it might have important effects on the structural response. As expected and highlighted by the HHs in Fig. 19.b he energy is concentrated around the first loading event with energy spread on a large range of frequencies. The presence of a force component coherent with the lower frequency identified within the IRF (i.e. 13 Hz), might signify that, despite the removal of the cylinder vibrations from the recorded signal, the original loading process was affected by the non-fully rigid nature of the experimental set-up. Hence the exerted hydrodynamic loading might not be exactly the same that would have occurred in a situation with a fully rigid structure.
Although the proposed method allows the identification of most of the main features of the wave loading, the description of the application point, calculated by the crude ratio between the overturning moment and the horizontal force, is still affected by some inaccuracies that are reflected in the gaps within the time series. The main reason is the level of noise that is present in the MAP solution. Indeed, the gaps in the application point time series (black line and grey area in Fig. 19c) correspond to the lower values of the identified wave force, so that the division between the moment and small force values provides unrealistic results. Furthermore, few values of the application point fall above the detected runup. This is associated with the inaccuracy of the image processing based measurement of the runup that sometimes fails to properly detect the high turbulent or aerated water mass. In Fig. 19.c, for the sake of visual rendering, all the values of the application point related to forces smaller than 2 N have been removed. Both Fig. 18b, d and Fig. 19a,c shows the credible interval around the MAP solution for the force, the moment and the application point; however, due to the small signal-to-noise ratio the posterior distributions are quite narrow, so that the shaded grey area is slightly obscured. Two examples of posterior distributions are presented in Fig. 19a,c for the instant related to the maximum identified force (i.e. 3.153 s) and are overlapped with the main time series within the smaller plot boxes at 3.153 s. Note how the uncertainty in the application point is quite large for small values of force, whereas it becomes relatively small for the main impulsive load, as a result of the stronger information carried by the data in the Bayesian process. Finally, Fig. 19.c also shows the average least-square solution for the application point as a black horizontal dotted line with the associated uncertainty, again quite small and barely visible in the figure. From the comparison between the time varying application point and the overall least-square solution it is clear the potential of the proposed analysis method and the need for a proper post-processing procedure for A. Antonini et al. the laboratory wave force time series.

Field data: Wolf Rock lighthouse
In order to show the capability of the proposed approach, the same methodology is applied to the field acceleration measurements recorded during Hurricane Ophelia (October 2017) on Wolf Rock lighthouse. The idea is to use the lighthouse as full-scale force transducer and reconstruct the wave load that induced the shaking of the structure. However, the field nature of this application requires more detailed characterisation of the dynamic behaviour of the structure, so that the modal parameters must be identified by means of a dedicated field modal campaign. The same approach can be applied to more common marine structures such as: vertical wall breakwaters, crown-walls, offshore wind turbines and offshore platforms.

Wolf Rock lighthouse: location, features and modal analysis
Wolf Rock lighthouse (49 • 56.72 ′ N -05 • 48.50 ′ W), Fig. 20 left panel, is situated 13 km offshore the most south-westerly point of UK, halfway between the Isles of Scilly and Land's End. It is one of the most exposed lighthouses in the British Isles, being surrounded by more than 35 m water depth on all sides but the south-east, Brownjohn et al. (2018); Raby et al. (2019a). The tower is composed of 70 granite courses, and extends to a height of 41 m from foundation to highest course. If the extent of the helideck is also considered (the first one constructed on top of a lighthouse, in 1973) it reaches a height of 43.1 m. Each granite course of the tower is subdivided into 16 sectors, each masonry course and sector being connected with their neighbours through vertical key and dovetail joints. The outside diameter reduces from a maximum of 12.68 m at the complete 2nd course to a minimum 5.18 m at the 68th course. The total volume of the granite is 1260 m 3 having a mass of 3350 t. The lower landing platform extends north-east for about 25 m and is covered by granite blocks about 0.15 m thick. More detailed descriptions of the Wolf Rock lighthouse can be found in Raby et al. (2019b) and Brownjohn et al. (2018).
A field campaign aimed to identify the lighthouse modal parameters such as modal masses, natural frequencies, damping ratios and mode shapes was performed in 2016 as part of the STORMLAMP project activities (Brownjohn et al. (2018)). During the two day campaign, both ambient and forced vibrations were recorded at the 8 + 1 floors (masonry tower and helideck) of the lighthouse. Orthogonal pairs of Honeywell QA-750 quartz-flex accelerometers were arranged at the inner wall of the masonry tower at the same compass bearing with respect to the lighthouse vertical axis, Fig. 20 right panel, while a shaker was located at the battery room and acted along both the x and y directions. The ambient vibration data were post-processed with standard Eigensystem Realisation Algorithm (ERA) (James et al., 1993) allowing the identification of natural frequencies, mode shapes and damping ratios. The forced vibrations were analysed with both Global Rational Fraction Polynomial (GRFP) (Richardson and Formenti, 1985) and circle fit (CFIT) functions (Kennedy and PANCU, 1947), additionally allowing the identification of the modal masses (Fig. 21).
Modal masses are most important for relating wave loading to response, but often they can be misunderstood because of the way they are linked to mode shape scaling. Therefore we define the modal mass as the integral with respect to height of mass weighted by squared horizontal modal ordinate (Brownjohn and Pavic, 2007), while the scaling sets the mode shapes to have unitary value at the level where the shaker is located (i.e. battery room). Accordingly, the physical response at this location is obtained by considering each mode as a SDoF system with this "unity scaled" value of modal mass. Fig. 21 shows the obtained results from the GRFP in which the first two identified mode shapes are presented. The mode shapes with large helideck ordinate have much larger modal mass. This is because the contribution to the modal mass calculation goes with the square of the modal ordinate. Furthermore, since the open helideck structure is practically transparent to horizontal loads due to breaking wave impacts, response of the masonry towers in these modes is expected to be relatively low.

Fig. 24.
Experienced inertial force at the battery room during one wave impact during the Hurricane Ophelia.

Hurricane Ophelia
Hurricane Ophelia, the tenth and final consecutive hurricane of the very active 2017 Atlantic hurricane season, was the strongest storm that affected the south-western UK and Irish coasts. Formed on 3rd October from a broad low-pressure area offshore the Azores, it began to strike the British and Irish coasts at the beginning of 12th October (Guisado--Pintado and Jackson, 2018). At the Wolf Rock nearest deep water node (49 • 56 ′ 4.7 ′′ N -5 • 49 ′ 46 ′′ W), available within the NORGASUG model (Boudière et al., 2013), the hurricane reached its maximum intensities in term of significant wave height on 16th October during the afternoon between 13:00 and 16:00 (H m0MAX = 7.15 m; T P = 13.3s, T S = 12.4s, T m = 11.2s and d P = 222 • N) and 21st October, during the morning, between 10:00 and 13:00 (H m0MAX = 7.55 m; T P = 14.7s, T S = 13.7s, T m = 12.4s and d P = 242 • N), Fig. 22. The following analysis will focus on the 16th October peak.

Remote acceleration acquisition system
During the same period a remote logging system, aimed to acquire the wave-induced acceleration, was installed in the lighthouse battery room, i.e. 7th floor. The system comprises a single JA-70SA triaxial servo accelerometer, therefore two horizontal (i.e. CH 1 and CH 2 ) and one vertical signals are available for the analysis. The vertical acceleration is negligible for the aim of the proposed analysis. CH 1 points 282 • N and CH 2 is perpendicular, i.e. it points 192 • N, Fig. 23 upper panel. The recorded accelerations are shown in Fig. 23 lower panels at different time scales, i.e. the entire day, the 6 h of the storm and the selected acceleration time series together with their HHs.
The original recorded CH 2 acceleration is rotated by approximately 30 • from the hindcast wave direction; thus, to extract the acceleration of the lighthouse only along the wave direction the simultaneous records of both perpendicular channels are iteratively rotated (by step of 1 • ) from 0 to 180 • clockwise. Moreover, the rotation analysis also allows the assessment of coherence between the acceleration measurements and the hindcast wave direction. The integral of the energy spectrum for both rotated signals is calculated and used as a proxy for the estimation of the impact direction. The maximum value of CH 2 energy spectrum integral is obtained for a clockwise rotation equal to 34 • , hence it can be argued that the wave generating the shaking of Wolf Rock lighthouse at 14:28 on 16th October was coming from 226 • N in close agreement with the hindcast wave direction. Thus in the following study the 34 • rotated CH 2 signal is considered, i.e. the considered alignment is 226 • N.
To identify the wave force generating the shaking at the battery room, the experienced inertia force should be described through the modal masses presented in Fig. 21 Brownjohn et al., 2018) and then used within the inverse process. In other words, we want to use the measured acceleration of an elastically linked mass to evaluate the force acting on the mass itself. In the lower panels of Fig. 23 the HHs for both CH 1 and CH 2 show that the impulsive wave load response is mainly concentrated within the second natural mode, i.e. around 6.8 Hz. Thus, under the assumption of linear behaviour of the structure, the experienced inertial force at the battery room can be calculated as the product between the acceleration and the modal mass related to the second natural mode (i.e. 436 t) as shown Fig. 24. To summarise this preliminary analysis of the acquired signal, impulsive wave load will be reconstructed by the application of the previously described inverse method considering the inertial force and the IRF related to the second natural mode.
The process follows the same steps adopted for the laboratory data analysis; however, since it is impossible to test the lighthouse with an impact hammer, the IRF for the second natural mode is reconstructed from its theoretical expression. The process requires the deconvolution between two homogenous signals, hence the IRF expressed in term of displacement is converted in terms of force by means of the product between the second mode modal mass and the second time derivative of the displacement IRF as in eq. (3). In this case the modal parameters identified through the field modal analysis are adopted within eq. (3). In this regard, it should be mentioned that we are assuming that the dynamic parameters identified through the dry modal analysis, i.e. when no water was in contact with the lighthouse, remain valid also during the wave impact; whereas in Brownjohn et al. (2019) the results of the non-stationary modal analysis during Hurricane Ophelia show the increasing of the damping ratio together with the decreasing of the natural frequencies during the period of maximum wave agitation. In such cases, it is reasonable to argue that the water surrounding the structure during the impact exerts additional damping and contribution to the inertia of the structure through the added mass. However, the limited number of observations and the uncertainty on the added mass do not allow the systematic use of these non-stationary parameters within the inverse process.

Field (almost) uninformative prior distributions
The lack of information about the individual wave that might have generated the analysed lighthouse shaking imposes the use of uninformative prior distributions about the intensity of the impulsive wave loading. Therefore, time constant normal prior distributions characterised by mean value equal to zero and relatively large standard deviation (i.e. 4 times the maximum calculated inertial force at the battery room) are adopted. Such a large value of the standard deviation makes the prior similar to a uniform distribution for the range of investigated values. In contrast, the time duration of the impulsive loading can be estimated with reasonable accuracy. We assume that the linear phase celerity is valid for the estimation of the velocity of the water mass that hits the lighthouse, hence, knowing the water depth (i.e. 35 m) and assuming that the largest waves are associated with the significant period equal to 11.75 s (Goda (2000)), the wave phase celerity is equal to 15.40 m/s. Moreover, applying Goda et al. (1966) impulsive loading duration model and considering the diameter equal to 12.68 m (the maximum at the base of the lighthouse) it is possible to make a preliminary estimate of the duration of the loading impulsive component of 0.41 s, hence the prior correlation function falls off with the same time scale. A detailed sensitivity analysis on the effect of the prior distributions on the final result is presented in appendix B for sake of completeness.

Field results
The result of the inverse process applied to the field data is presented in Fig. 25, where the overall inertia force due to the 16th October event is shown as a dotted light blue line. The reconstructed force experienced at the battery level due to the wave impact is presented with a solid black line (the credible interval is not visible due to the scale) the lower panel shows its HHs. The lower left box shows the detail of the posterior distribution related to the maximum reconstructed force value, while the upper right box is the enlargement of the reconstructed peak force.
Not only for the well-controlled laboratory data, but also for the more complex field data it can be said that the inverse method works properly and a large amount of the structural dynamic effects are removed from the signal, allowing a clear description of the wave force features. Weak background oscillations remain in the final reconstructed force due to the theoretical nature of the adopted IRF that does not provide an ideal kernel for the deconvolution of complex field data. Some of the dynamic features of the lighthouse are not perfectly removed from the final solution, as shown by the oscillations characterised by a frequency slightly larger than the first natural mode between 14:28:40 and 14:28:42. However, despite this spurious less energetic component in the final result, the impulsive components are properly captured.
Similar features that were previously identified in the laboratory result are also detectable for the field outcomes. The two close peaks characterising the first impulsive loading component are properly reconstructed and likely to be due to the impact of the first front and the following deceleration of the falling water mass. Although the initial impulsive component does not show the same characteristics of a fully breaking wave loading, its slightly longer character remains rather similar to a slamming force. The energy is concentrated within a few tenths of second, between the beginning of the rise time and the end of the second peak. The instantaneous frequency shows the typical feature of the impulsive load, with energy content spread over a rather large frequency band, ranging between 0.8 and 40-50 Hz. While the laboratory results showed clearer concentration of energy at the beginning of the rise time, here the maximum energy is detected at the beginning of the second peak; however, a sudden increase of the instantaneous frequency is also visible at the beginning of the overall impulsive load, highlighting again the importance of this initial part of the load.
The first, double peak, impulsive component lasts for about 0.06 s (slightly before 14:28:39) that in turn makes the nature of the load nearly dynamic with respect to the behaviour of the structure (i.e. the ratio between the impulsive loading length and the natural period is larger than 0.25 (Chen et al., 2019;Oumeraci and Kortenhaus, 1994)). This aspect is particularly relevant for the accurate modelling of the response of the structure under wave action. Indeed, if the first mode would have been the most significant one, the structural analysis could have been carried out according to the impulse theory presented by Chen et al. (2019), whereas, in this situation, the time-varying nature of the impulsive wave load must be taken into account. Proceeding along the reconstructed incident wave loading, a third less intense peak is present and highlighted by the sudden increase of the instantaneous frequency at about 14:28:39.5. The similarity between this third peak and the identified one in the laboratory result is clear. However, a slightly different time sequence, i.e. time lag between the main impulsive component and the third peak, is also evident and likely to be due to the uncertain condition at the base of the lighthouse and the induced effects on the wave breaking process.
Up to this point in the analysis, the adopted methodology has virtually assumed the application of the incident wave force at the battery level, i.e. 32.6 m above the first full course of the lighthouse or 34.7 m above the chart datum. Evidently, the wave action could not be directly applied to such a high level. One last step is therefore required to translate the application point from the battery room to a more realistic lower application point, and accordingly re-scale the wave force intensity. This final step can be achieved by knowing the modal ordinates at different elevations of the structure. Within the proposed analysis, the battery room modal ordinate was always kept equal to 1. Hence, in order to reconstruct the force at a lower level, the reconstructed force time series should be divided by the modal ordinate (in this case related to the second mode) corresponding to the estimated application point. However, since no direct measurements of the incident individual wave height are available for Wolf Rock lighthouse, a preliminary estimation of the wave force application point is carried out through the available methods in the literature. The statistical distribution proposed by Battjes and Groenendijk (2000) is adopted to describe the individual wave heights at the toe of the lighthouse (Fig. 26  upper panel), while the application point is calculated under the assumption of uniform vertical pressure distribution exerted along the upper 46% (i.e. curling factor equal to 0.46) of the asymmetric wave crest described by means of Hansen (1990) method ( Fig. 26 middle  panel).
Once the application point is known, the modal ordinate value is also known and can be used as the scale factor for the intensity of the reconstructed force (Fig. 26 lower panel). From Fig. 23 the analysed impact is the 2nd highest in 3 h (i.e. 13:00 to 16:00) characterised by almost uniform H S around 7.15 m. Hence, considering the mean period (T m ) associated with the underlying Jonswap spectrum, 812 events can be estimated and accordingly an exceedance probability of 0.21% for the second-highest event can be calculated. Fig. 27 shows the final results of the overall inverse force reconstruction process. Due to the lack of direct measurement of the incident waves, the results about which wave had generated the recorded shaking are still affected by some uncertainties, however, they lay the foundation for a process-based assessment of the structure, where no empirical formulae are required to describe the incident wave loading.

Discussion and conclusion
This work intends to make progress in the application of the inverse method to reconstruct wave forces exerted on marine structures, providing a sound framework for a large number of field and laboratory applications. The presented methodology is based on linear theory and therefore assumes the elastic behaviour of the investigated structure, an aspect that should first be checked when the method is applied to field measurements under extreme wave actions. The paper provides a comprehensive presentation of the method by addressing three different applications: a theoretical one, where the impulsive load is reconstructed from the response of a theoretical single degree of freedom system; a laboratory application, where the force exerted by a broken wave is identified from the force measured on a vertical cylinder upon a shoal; and a field application, where a wave force exerted on Wolf Rock lighthouse is described from the accelerations measured during the Hurricane Ophelia. In the following, the main aspects of the three applications are discussed and some conclusions gathered. The theoretical example shows, in a simplified and controlled case, the main issues related to the inverse process. The effect due to the (inevitable) presence of noise within the real signal is highlighted and analysed in detail. The methodology is applied and the results compared with the Tikhonov regularization and the truncated Fourier transform. Although the truncated Fourier transform is rather time-efficient (the computational time is way smaller than the other methods) and effective to reconstruct the unknown force, the method requires the subjective selection of the threshold to overcome the sensitivity of the results to small changes in the input records (noise) and the results show loss of resolution and some spurious oscillations. The Tikhonov method is effective for the reconstruction of the essential features of the impulsive input signal, the truncated singular values show their influence through the introduction of spurious oscillations and loss of resolution in the final solution. On the other hand, the Bayesian methodology allows a detailed reconstruction of the input signal. The improvements due to the introduction of process-based informative prior distributions are evident on both the final solution and the narrow posterior distribution; indeed the sharp nature of the impulsive wave loading is reconstructed without the presence of substantial spurious oscillations. The theoretical example is also available as additional material in the form of MATLAB® procedures.
For the laboratory data application, the use of an impact hammer aimed to define the experimental force IRF is presented. The related illconditioning issues are tackled through the combined use of the leastsquares method and QR decomposition. The identified experimental IRF is then applied to five regular waves. The main features of the load due to the broken waves are captured in detail by the method, allowing an accurate time-varying description. The load is characterised by an initial impulsive component constituted by two consecutive peaks and a delayed one. The first peak is mainly due to the impact of the first aerated front, while the second is due to the sudden deceleration of the falling water mass at the base of the structures. The third peak has a smoother shape, is less intense and is longer than the first ones. It is due to the impact of the remaining water mass carried by the individual wave. Both wave force and overturning moment are reconstructed with a good level of accuracy, allowing the identification of the point of application of the force. The application point time series is partially affected by the remaining small spurious oscillations in the final solution. The effect is particularly evident for low force intensity. In this condition, the resulting time series is overcome by the spurious oscillations making part of the application point results unreliable. However, the application points for the impulsive loading conditions are properly captured. The goodness of the final result is further corroborated by the comparison between the force, overturning moment and force application point with the detected cylinder runup. Although the overall analysis provides trustful results some uncertainties remain on the effects of the structure compliance and its movements during the interaction with the water mass, however, they are hardly quantifiable and assumed to be negligible.
The application of the Bayesian method to the field vibration data is slightly more complex. An accurate dynamic characterisation of the structure is required in order to identify the main dynamic parameters and mode shapes. The data is initially pre-processed to indirectly identify the direction of the incident wave generating the recorded shaking. The identified wave direction agrees with the direction coming from the hindcast model, validating the adopted analysis procedure. For Wolf Rock lighthouse only the second identified mode seems to respond to the impulsive wave loading. Based on this finding the deconvolution process is based on the theoretical IRF of the same mode. The final result allows the identification of the time-varying nature of the wave load. The same features identified for the laboratory results are also detected for the field data. Three peaks characterize the reconstructed wave loading and can be argued that they are generated by same physical phenomena observed in the laboratory experiments. The field data do not allow the identification of the overturning moment or the direct measurement of the incident wave; accordingly also the wave force application point is not directly described. To overcome this lack of information a statistical description of the possible incident waves is performed in order to estimate the constant application point of the force and then rescale the intensity of the reconstructed force via the inverse of the modal ordinate. According to the second natural mode, the lower the wave force application point, the larger should have been the force intensity. From the recorded vibrations a preliminary estimation of the probability of exceedance of the analysed event is inferred equally to 0.21%, accordingly the generating individual wave and force application point are calculated.
Overall, the proposed methodology allows the reconstruction of the wave force directly from structural dynamic measurements, laying the foundation to analyse unclear physical phenomena such as breaking and broken wave loading on rigid structures. The method can be extended to multiple degrees of freedom as well as to structures that respond to the wave loading with a combination of multiple natural modes. Here, it has been applied to a laboratory cylinder or to an offshore rock lighthouse; however, we expect that the same procedure should be effective also for more common marine structures such as caissons, crown-walls, wind turbine monopile and offshore platforms. For future field applications, we strongly advise the planning of a measurement campaign where the simultaneous record of the structural vibration and individual incident wave heights are considered.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Laboratory unit-impulse response function
A theoretical unit-impulse response function (IRF) can be calculated as in eq. (17):

FFT(System output) FFT(System impulsive input)
) ∫ t1 0 System impulsive input dt (17) where t 1 is the length of the available record, FFT is the fast Fourier Transform and IFFT is its inverse. Moreover, it must be noted that the ratio in the argument of the IFFT is just the time domain deconvolution between the system output and input that can be resolved by the matrices operation presented in eq. (16). In the proposed laboratory application the overall signal presented in Fig. 15 is divided in three parts of the same length equal to 2.06 s that are later used to calculate three IRFs that in turn are averaged to obtain the final IRF used in the inverse method.
In the specific case of the laboratory force measurements the IRF is calculate as in eq. (18): where F output (t) is the force time series as recorded by the force transducer and F input (t) is the applied force as recorded by the impact hammer. If the same is applied to the calculation of the moment IRF the product of the force and its application point should be considered as in eq. (19): where b (t) is the time series of the application point of the applied force via the impact hammer, but also the application point comprising the moment measured by the moment transducer. b (t) is the same in both numerator and denominator of the IFFT argument, so the value can be cancelled out and therefore considered negligible for the calculation of the moment IRF. Obviously, the normalising integral should be redefined as in eq. (18), however, its result is a scalar that is only needed to normalise the IRF and therefore it only affects the amplitude of the IRF oscillations whereas the other features remain unchanged. and eq. (19) respectively. The assumed constant application point is equal to 0.525 m from the origin of the reference system shown in Fig. 11. Fig. 28 shows the comparison of the two IRFs calculated by explicitly taking into account both impact hammer test response for the force and the moment under the assumption of a distance between the reference system origin and the impact hammer application point equal to 0.525 m. This distance was not measured during the impact hammer tests. Indeed, despite the fact that we roughly applied the impact hammer around the area affected by the wave action, we could not precisely measure the position and had to estimate it afterwards, introducing an additional uncertainty in the model. Both upper and lower panels show a high-frequency component for the force IRF that is not present for the moment IRF, approximately around 500 Hz as highlighted in the spectrum (lower panel). This high-frequency component can only be due to the internal transducer set-up, as both the force and moment have been measured by the same integrated transducer. In light of this difference in the frequency contents and the uncertainty due to an erroneous impact hammer application point, we had to make a choice about what could be considered negligible. The choice was to proceed only with the force IRF for better control we could have on the final solution and the reasonable theoretical background provided above (Fig. 29).

Appendix B. Prior distribution sensitivity analysis
Due to the limited prior process-based knowledge on the wave event generating the shaking of Wolf Rock lighthouse on October 16, 2017 almost uninformative uniform prior distributions have been applied for the proposed analysis. However, to investigate the effect of the different prior distributions parameters, a sensitivity analysis on the final results has been performed. Being the impulsive load component, the most interesting part from a structural point of view, we focused the sensitivity analysis on three features describing this component, i.e. the maximum force, the time duration and its impulse quantity. The wave slamming component is defined as the reconstructed wave load time series between the zero-down and -up crossing points nearest to the maximum reconstructed force, while its underlying impulse is defined by the integral of the force-time series between these two points, Fig. 29.   Fig. 29. Example of Wolf Rock lighthouse reconstructed wave force and the adopted parameters for the prior distributions sensitivity analysis The sensitivity analysis focused on two main parameters. The standard deviation (std) for the prior normal distribution, expressed in term of percentage of the maximum absolute value of the recorded inertial force (Fig. 24), with values between 10% and 600%. While the second parameter is the correlation time lag, that has been varied from 0.026 s to a maximum of 0.78 s. The results for both parameters are presented in Fig. 30, where the right column shows the results related to the prior distribution std and the left column the results related to the correlation time lag. It is evident that, except for small values of the std the selection of different parameters do not largely affect the final result. Indeed, for std values larger than 2 times the maximum recorded inertial force no evident variation in the final result can be detected, while also the correlation time lag is not a critical parameter being quite weak the trend within the results. Therefore, we adopted middle values for both the parameters, keeping the prior std equal to 4 times the maximum recorded inertial force and the correlation time lag around 0.41s also according to the preliminary estimation of the impulsive component.