# An indirect data-driven model-updating framework to estimate soil-pile interaction parameters using output-only data

Andreas Ioakim, Szymon Greś, Michael Döhler, Luke J Prendergast, Eleni Chatzi

<figure><img src="image_1.png" /></figure>

To cite this version:

Andreas Ioakim, Szymon Greś, Michael Döhler, Luke J Prendergast, Eleni Chatzi. An indirect data-driven model-updating framework to estimate soil-pile interaction parameters using output-only data. Engineering Structures, 2025, 328, pp.119699. 10.1016/j.engstruct.2025.119699. hal-04902544

HAL Id: hal-04902544

https://inria.hal.science/hal-04902544v1

Submitted on 21 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An indirect data-driven model-updating framework to estimate soil-pile interaction parameters using output-only data

Andreas Ioakim<sup>a,*</sup>, Szymon Gres<sup>b,c</sup>, Michael Döhler<sup>d</sup>, Luke J. Prendergast<sup>a</sup>, Eleni Chatzi<sup>b</sup>

<sup>a</sup>Department of Civil Engineering, University of Nottingham, Nottingham, NG7 2RD, United Kingdom

<sup>b</sup>Department of Civil, Environmental and Geomatic Engineering, ETH Zürich, 8093 Zürich, Switzerland

<sup>c</sup>Department of Civil and Architectural Engineering, Aarhus University, 8000 Aarhus, Denmark

<sup>d</sup>Univ. Gustave Eiffel, Inria, COSYS/SII, I4S, Campus de Beaulieu, 35042 Rennes, France

## Abstract

This paper presents a data-driven model updating framework to estimate the operational parameters of a laterally-impacted pile. The goal is to facilitate the estimation of soil-pile interaction parameters such as the mobilized mass and stiffness, as well as geometrical data such as embedded pile length, using output-only information. Accurate knowledge of mass, stiffness, and pile embedded length is essential for understanding foundation behavior when developing digital-twin models of structures for the purpose of damage detection. The method first employs subspace identification to determine modal parameters and quantifies their uncertainties using output-only data. The covariance matrix adaptation evolution strategy (CMA-ES), a stochastic evolutionary algorithm, is subsequently used to update the model. The effectiveness of the approach is demonstrated through its application to numerical models in this paper, to quantify errors, and subsequently to data from a documented full-scale field test of a pile subjected to an impact load. The work underscores the potential of statistical updating in advancing the accuracy and reliability of soil-structure interaction parameter estimation for systems where only output data might exist.

**Keywords:** Model Updating; Structural Health Monitoring; Soil Structure Interaction; Winkler Model; Length Estimation;

## 1. Introduction

Analysing the response of piled structures commonly makes use of simplified models characterising soil-pile interaction, such as Winkler models, which model soil-structure behaviour using a series of distributed springs and masses attached to beam elements [1–5]. While the Winkler simplification is reasonably widely-adopted, the real soil-structure interaction behaviour is a complex three-dimensional phenomenon, and the dynamic response of soil depends on a variety of factors such as particle size, geological history, age of the deposit, degree of cementation, and the effective weight and stiffness of the soil surrounding the pile, among

<footnote>*Corresponding author; E-mail address: Andreas.Ioakim@nottingham.ac.uk</footnote>
other factors [6]. Simplified Winkler-type models offer an estimate into the behaviour of the real system, but if accurate behaviour is required, simplified models can be improved using model-updating approaches [7, 8].

Structural Health Monitoring (SHM) refers to monitoring the integrity of structural systems using mea-
surements collected from structures in-operation [9]. SHM covers an extensive range of topics, including
data-driven modal analysis, structural damage diagnosis, prediction of remaining useful life, model-updating
for parameter estimation, among many others [10]. For inverse problems, where the outputs are measured
but the system properties are unknown, data-driven modal identification can be used to estimate the modal
parameters of the monitored system [11, 12]. Model-updating methods can subsequently be used, whereby
the properties of the reconstructed model are calibrated based on the information obtained from the modal
identification [13], allowing for improved accuracy of the model compared to the reference model.

The primary objective of modal identification is to estimate the modal parameters (eigenfrequencies,
damping ratios, and mode shapes), which characterize the dynamic behaviour of the monitored physical
system. Depending on the availability of measured forces exciting the system, output-only or input/output
methods are used for the identification, where a wide range of well-established methods in time or in
frequency domain are available [11, 14, 15]. In particular, subspace-based system identification methods
are known for their ability to handle noisy data obtained from impact dynamic testing [16–18], statistical
consistency [19], and the ability to provide confidence bounds on quantities of interest [20–24].

Calibrated models are invaluable for understanding system behaviour, or for the purpose of monitoring
the structural state. Model-updating can be categorized into two main groups of methods: direct methods,
which employ motion equations and mode orthogonality properties to derive an analytical solution; and
indirect methods, which adjust model parameters based on experimental data through repeated calculations
[25]. Direct methods, praised for their computational efficiency and avoidance of divergence issues, aim for
an exact match with the measured data. However, they require precise, noise-free measurements and com-
plete modal data, which may not always be available for large-scale structures and are generally impractical
to achieve [26, 27]. Moreover, the obtained updated matrices might not always maintain physical rele-
vance [28, 29], i.e. the responses might match but the updated system might not be physically-admissible.
Given these limitations, indirect methods have been developed to address the challenges inherent to direct
methods [25]. Indirect methods involve calculating and minimizing an error function repeatedly to up-
date model parameters, and come in various forms, including sensitivity-based methods [30, 31], Frequency
Response Function-based methods [32–34], Bayesian theorem-based methods [35–37], and various compu-
tational techniques [38–40]. These approaches offer more flexibility and can accommodate the complexities
and imperfections of real-world data and large-scale system analysis, providing a more practical route for
model-updating in structural engineering.

This paper presents an application of an indirect model-updating scheme from Gres et al. [39] to
determine the embedded length of foundation piles, along with the mass of mobilized soil and the stiffness under ‘unknown’ input loading. Estimating the embedded length of piles has been explored in the past by Refs.[41–46], mostly using guided wave theory. The employed methodology in this paper builds on the study of Ioakim and Prendergast [34], which used input-output data. In the present work, output-only data is used and the method focuses on combining the data and model information with a stochastic optimization of the sought-after model-based parameters. The objective is to demonstrate the capability of this optimization approach in accurately predicting the distribution of essential soil-pile interaction parameters, without the need to measure the input force time-history. The validity and applicability of the study are evaluated through both numerical simulations and experimental applications. The proposed approach can serve as a meaningful tool for offshore wind turbine foundation analysis, for example, whereby knowledge or direct measurement of the input loads is typically not possible.

The structure of the paper is as follows: Section 2 provides an overview of the model updating framework. This section details the monopile system modeling using 1-D modelling, followed by the formulation of the objective function for model updating and the CMA-ES optimization strategy. Section 3 presents a numerical application to test the framework. Section 4 evaluates the model updating procedure using data from a documented large-scale field test on a pile foundation. Final remarks on the effectiveness of the framework are discussed in Section 5.

## 2. Methodology

The objective of this work is to develop an approach that can estimate the physical properties of both the soil and the pile in a pile-soil interaction problem, where the excitation applied is assumed unknown. The estimates of physical parameters are obtained through stochastic optimization of an objective function that confronts the modal parameter estimates of the pile within their estimation uncertainties to the corresponding parameters of the pile model. This is based on the evolutionary optimization strategy devised in Ref. [39].

Hereafter the parametric model of the pile foundation is described and the optimization procedure is outlined.

### 2.1. Pile-soil interaction modelling

The dynamics of the monitored parametric linear time-invariant system with $d$ degrees of freedom (DOF) are described by the differential equation of motion. The governing differential equation (1) is linear in nature, as the system’s response depends linearly on the external forces and state variables. The assumption of time-invariance is appropriate because the system matrices do not explicitly change with time during the analysis [1, 3].

$$M^{\theta}\ddot{q}(t) + D^{\theta}\dot{q}(t) + K^{\theta}q(t) = u(t) \quad (1)$$
where $t$ denotes continuous time, $M^{\theta}$, $D^{\theta}$, and $K^{\theta} \in \mathbb{R}^{d \times d}$ denote the mass, damping, and stiffness matrices, respectively and $\theta$ is the system parameter. For the considered problem, the system parameter $\theta$ holds the embedded pile length $L_1$, which is the geometrical data of interest, and the mass and stiffness weightings ($w_m$) and ($w_k$) that scale the soil's mass and stiffness profiles, respectively. These weightings are used to define the soil's mobilized mass and stiffness profiles [6]. The vectors $q(t)$ and $u(t) \in \mathbb{R}^d$ denote the continuous-time displacements and the unknown external forces, respectively.

The pile-soil system is modelled using a linear Winkler model, which utilizes 6 DOF Euler-Bernoulli beam/pile elements to characterize the pile, paired with 2-DOF mass-spring elements to represent the soil behavior. The model is capable of simulating the behaviour of soil-pile interaction at low strains, whereby the response of the system can be assumed linear [47, 48]. This is appropriate to systems oscillating with low-amplitude vibrations. The local mass and stiffness matrices for the pile elements can be calculated by simplifying the Timoshenko beam theory by ignoring the effects of shear deformation and rotatory inertia [49]. The 2-DOF soil springs are connected to the lateral DOF of each node within the embedded section of the pile. The pile's base is fixed in the vertical direction to exclude axial deformations from the system. This boundary condition effectively eliminates axial DOF contributions, while preserving the intended lateral behavior of the system. Each pile element of the pile is associated with two corresponding mass-spring elements that simulates the lateral mass and stiffness provided by the soil, extending up to the $N_s$th spring, which represents the topmost spring element. Specifically, the pile is discretized into a finite number of pile elements, and each $m$th pile element has an associated $n$th and $(n+1)$th mass-spring element. A schematic of the system is illustrated in Figure 1, which illustrates the physical system and the numerical spring-beam model representation.
<figure><img src="image_3.png" /><figcaption>Figure 1: Schematic of the pile model. (a) Physical model, (b) Numerical representation.</figcaption></figure>

The pile is assumed to be embedded in a single layer of cohesionless sand. The model captures the pile's flexural rigidity by incorporating the relationships between the subgrade reaction modulus and the pile material's elastic properties after Ref. [50].

$$K_S(z) = \frac{0.65 E_S(z)}{D(1 - \nu_S^2)} \sqrt{\frac{12(E_S(z)D^4)}{EI}}, \quad (2)$$

where, $E_S(z)$ represents the soil's small-strain Young's modulus (N m⁻²) as a function of the depth $z$, $\nu_S$ is the soil's Poisson's ratio, $D$ is the pile's width (m), and $E$ and $I$ denote the pile's Young's modulus (N m⁻²) and moment of inertia (m⁴), respectively. The small-strain Young's modulus of the soil is calculated as

$$E_S(z) = 2G_S(z)(1 + \nu_S), \quad (3)$$

where $G_S(z)$ is the shear modulus of the soil (N m⁻²) as a function of $z$. Given that $K_S(z)$ has units of N m⁻³ and the springs provide one-dimensional stiffness (N m⁻¹), the stiffness for each nth spring element is determined as

$$k_n = K_S(z)DL_m, \quad (4)$$

where $L_m$ is the tributary length of the associated mth pile element, and $K_S(z)$ is defined in Eq. (2). In the embedded layer, the pile nodes are attached to lateral springs, while those above ground remain free. To
simulate the dynamic influence of the mobilized mass of soil, lumped masses calculated as a fraction of the pile mass are distributed along the soil springs. Masses are added along the upper quarter of the embedded length, which is based on the premise that deeper masses exert a negligible effect on the first mode shape of a laterally-oscillating pile, leading to their exclusion [6]. The mass for each nth spring is therefore a function of the pile's mass

$$M_n = \frac{M_{\text{Pile}}}{Q}, \qquad (5)$$

where $M_{\text{Pile}}$ is the pile's mass, and $Q$ represents the number of springs in the top quarter of the pile.

## 2.2. Design of the objective function

The true modal parameters of system described in Eq. (1) are obtained from the eigenvalues $\mu_i^\theta \in \mathbb{C}$ and eigenvectors $\Psi_i^\theta \in \mathbb{C}^d$, which are the solution of the eigenvalue problem

$$((\mu_i^\theta)^2 \mathcal{M}^\theta + \mu_i^\theta \mathcal{D}^\theta + \mathcal{K}^\theta) \Psi_i^\theta = 0, \qquad (6)$$

and yield the natural frequencies $f_i^\theta = \frac{|\mu_i^\theta|}{2\pi}$ and the mode shapes $\psi_i^\theta = L\Psi_i^\theta \in \mathbb{C}^r$ at the observed DOF for a given sensor type, where $L \in \mathbb{C}^{r \times d}$ maps the eigenvectors to the mode shapes that can be obtained from measurements at r sensors, respectively.

$\hat{f}_i$ and $\hat{\varphi}_i$ are denoted as the respective estimates of the i-th natural frequency and the i-th mode shape. The deployed system identification method is assumed consistent, i.e., $\hat{f}_i$ and $\hat{\varphi}_i$ converge to their true values $f_{i,*}$ and $\varphi_{i,*}$ when the data length $N$ goes to infinity. This assumption is leveraged by multiple system identification methods, for example stochastic subspace identification (SSI) methods [16, 21, 51], where these parameters can be identified, e.g., from acceleration, velocity, displacement, or strain measurements. There is an inherent assumption that the system is subject to ambient white noise excitation. It should be noted that depending on the physical case considered, this assumption might not be fully valid. For example, this assumption might not be exactly true in the case of operating offshore wind turbines supported on monopiles. The modal parameter estimates have been proved to be asymptotically Gaussian for SSI methods [22, 52], satisfying

$$\hat{f}_i \approx \mathcal{N}(f_{i,*}, \frac{1}{N}\sigma_{f_i}^2), \quad \text{and} \quad \begin{bmatrix} \Re(\hat{\varphi}_i) \\ \Im(\hat{\varphi}_i) \end{bmatrix} \approx \mathcal{N}\left( \begin{bmatrix} \Re(\varphi_{i,*}) \\ \Im(\varphi_{i,*}) \end{bmatrix}, \frac{1}{N}\Sigma_{\varphi_i} \right), \qquad (7)$$

where $\Re(\cdot)$ and $\Im(\cdot)$ express the real and imaginary parts of a complex variable, $\mathcal{N}(\mu, \Sigma)$ denotes a Gaussian distributed variable with mean $\mu$ and covariance $\Sigma$, and $\sigma_{f_i}^2 \in \mathbb{R}$ and $\Sigma_{\varphi_i} \in \mathbb{R}^{2r \times 2r}$ are the asymptotic covariance of $\hat{f}_i$ and $\hat{\varphi}_i$, respectively. The computation of the respective covariance estimates $\hat{\sigma}_{f_i}^2$ and $\hat{\Sigma}_{\varphi_i}$ from the same data set as the modal parameter estimates can be found, e.g., in Refs.[20–22]. Since $\hat{f}_i$ is approximately Gaussian distributed, its confidence interval is symmetric and can be centered around the computed estimate. As such, it will cover a range of plausible values that include the true value of the
parameter with some given confidence level $\gamma$, e.g., $\gamma = 0.9544$ for the $\pm 2\sigma$ bound. At this value, the corresponding confidence interval

$$c_{\hat{f}_i} = [\hat{f}_i^{\min}, \hat{f}_i^{\max}] \qquad (8)$$

encompasses the true value $f_{i,*}$, where $\hat{f}_i^{\min} = \hat{f}_i - 2 \cdot \frac{1}{\sqrt{N}} \hat{\sigma}_{f_i}$ and $\hat{f}_i^{\max} = \hat{f}_i + 2 \cdot \frac{1}{\sqrt{N}} \hat{\sigma}_{f_i}$ with 95.44% probability, which is hereafter denoted as the 95% confidence interval. If the natural frequency of the model $f_i^{\theta}$ is inside $c_{\hat{f}_i}$, then it can be inferred that both $\theta$ and $\theta_*$ are statistically close. Conversely, when $\theta$ is different from $\theta_*$, the corresponding model frequency $f_i^{\theta}$ is likely to lie outside the confidence interval of the estimated frequency.

An objective function for the model optimization is expressed by the difference between the modal parameter estimates and their model-based counterparts. It is hereafter assumed that the modal parameter pairing is accomplished with statistical criteria [53]. Denote $\Delta_{f_i}^{\theta}$ as the normalized difference between the estimated natural frequency $\hat{f}_i$ and the natural frequency $f_i^{\theta}$ from the model under parameter $\theta$, with

$$\Delta_{f_i}^{\theta} = \begin{cases} 1 - \frac{f_i^{\theta}}{\hat{f}_i}, & \text{if } \hat{f}_i^{\min} \le f_i^{\theta} \le \hat{f}_i^{\max} \\ 1, & \text{otherwise,} \end{cases} \qquad (9)$$

where a penalty of 1 is introduced when $f_i^{\theta}$ is outside the 95% confidence interval of $\hat{f}_i$. For the comparison of the mode shape estimate $\hat{\varphi}_i$ and model-based mode shape $\psi_i^{\theta}$, the Modal Assurance Criteria (MAC) is used. If $\theta = \theta_*$, the MAC between both mode shapes is close to 1 and tends to 1 when the data length $N \to \infty$, since $\hat{\varphi}_i$ is a consistent estimate. Consequently, if the MAC between both mode shapes tends to 1, the parameter vector $\theta$ is a valid candidate for $\theta_*$. Let $\Delta_{\text{MAC}_i}^{\theta}$ be the difference between 1 and the respective MAC, defined as

$$\Delta_{\text{MAC}_i}^{\theta} = 1 - \text{MAC}(\hat{\varphi}_i, \psi_i^{\theta}), \qquad (10)$$

where

$$\text{MAC}(\hat{\varphi}_i, \psi_i^{\theta}) = \frac{|\hat{\varphi}_i^H \psi_i^{\theta}|^2}{\hat{\varphi}_i^H \hat{\varphi}_i \psi_i^{\theta H} \psi_i^{\theta}}.$$

The objective function $F(\theta)$ is expressed as the sum of the respective differences $\Delta_{f_i}^{\theta}$ and $\Delta_{\text{MAC}_i}^{\theta}$ for all considered modes $i = 1, \dots, n_m$, as

$$F(\theta) = \sum_{i=1}^{n_m} \left| \tilde{\Delta}_{f_i}^{\theta} \right| + \sum_{i=1}^{n_m} \Delta_{\text{MAC}_i}^{\theta}. \qquad (11)$$

Based on the confidence intervals (8), the acceptance region

$$\tilde{\Theta} = \{\theta : f_i^{\theta} \in c_{\hat{f}_i} \text{ for all } i = 1, \dots, n_m\} \qquad (12)$$

can be defined, which comprises all parameter vectors $\theta$ that yield model frequencies within the confidence intervals of their estimated counterparts. Thus, $\tilde{\Theta}$ comprises the statistically acceptable solutions for $\theta_*$ with regards to the considered frequencies, and a stopping criterion of the optimization search can be formulated as $\theta \in \tilde{\Theta}$.
## 2.3. Model optimization strategy

The solution $\theta^{\text{sol}}$ is obtained via minimization of $F(\theta)$

$$ \theta^{\text{sol}} = \arg\min_{\theta} F(\theta), \qquad (13) $$

which is achieved with an application of the CMA-ES [54]. Starting with some initial value $\theta = \theta_{\text{init}}$, the CMA-ES algorithm consists in generating $\lambda$ model candidates $\theta_j^g$, $j = 1, \dots, \lambda$, in each population $g$, by sampling a multivariate Gaussian distribution. The sampling is carried out as

$$ (\theta_j)^{g+1} = m^g + \epsilon_j, \quad \epsilon_j \sim \sigma^g \mathcal{N}(0, C^g), \qquad (14) $$

where $j = 1, \dots, \lambda$ and $m^g$ is a weighted mean of the model candidates $(\theta_j)^g$ in the parent generation. Eq. (14) represents a mutation and recombination into offspring, for which the CMA-ES algorithm adapts the parameters $C^g$ and $\sigma^g$ in each generation. The covariance matrix $C^g$ denotes the covariance of the search distribution and $\sigma^g$ denotes the scaling factor adapted (minimized) in the optimization procedure. Consequently, the optimization continues and the best parent solutions replace the offspring until it converges to a solution. For CMA-ES, the covariance matrix $C^g$ is incrementally updated with rank-one matrices representing the direction between the best parent solutions at two consecutive generations, such that the likelihood of previously successful search steps is increased [54]. For the convergence to a solution, a stopping criterion is included in the algorithm that is adapted to the acceptance region (12), which is defined by the confidence intervals of the modal parameter estimates obtained from the sensor data. Once inside the acceptance region, there is no need to further minimize the objective function, avoiding unnecessary additional computations.

Figure 2 presents a flowchart of the model updating framework, demonstrating how the measurement data are linked with the SSI to obtain identified frequencies and mode shapes, coupled with their uncertainties. These identified modal properties are then combined with the computed modal properties from the numerical model to formulate the objective function used in the CMA-ES optimization.

<figure><img src="image_4.png" /><figcaption>Figure 2: Model updating flowchart.</figcaption></figure>
### 3. Numerical Application

In this section, the efficacy of the proposed framework in estimating soil-pile interaction parameters, including mobilized soil mass and stiffness, and embedded pile length, is evaluated. To undertake this demonstration, known values of soil mass and a profile of soil stiffness are altered by applying weightings $w_m$ and $w_k$, and assessing how well the method can estimate these weightings, i.e. converge on the actual soil mass and stiffness profiles (see Ref.[6] for more details on this approach). Weightings $w_m$ and $w_k$ are defined as scalars that adjust a given spring's mass and stiffness values. The optimization methodology is applied to a set of target pile models (with varying geometries and soil profiles) to assess its accuracy in parameter estimation. For more information on the development of models for this type of application, see Ioakim and Prendergast [34].

Three target models (pseudo-experimental piles subjected to simulated impact tests) are developed to represent various types of piles, which are referred to as target models throughout this study. The schematic of the pile configuration is shown in Figure 3. As mentioned earlier, the mobilized mass of the soil is assumed as concentrated in the top quarter of the pile, in keeping with the mobilisation expected under vibration in the first mode. For simplicity, only some nodes are displayed in the schematic, with intervals of 1 meter. Table 1 provides specific parameter values used to create these models. To examine applications to both flexible and rigid piles, two radii ($r$) with annular cross-sections and corresponding wall thicknesses ($t_0$) are explored: $r = 0.5$ m with $t_0 = 0.025$ m, and $r = 2$ m with $t_0 = 0.1$ m. Each model is constructed by discretizing the pile into finite elements of 0.1 m in length. For all target models, the total pile length ($L_0$) is defined as the sum of the embedded length plus an additional 1 m stick-up ($L_0 = L_1+1$ m). These parameters are chosen to represent a range of foundation geometries typically found in various applications, including the simulation of monopiles used in offshore wind turbine foundations [55]. Additionally, to simulate real field scenarios where sensor data might be influenced by environmental factors, noise is introduced into the acceleration data.
<figure><img src="image_5.png" /><figcaption>Figure 3: Pile configuration for numerical application.</figcaption></figure>

The soil's shear modulus profile is assumed as constant, analogous to an over-consolidated soil deposit, and is set at 50 MPa. Fig. 4 illustrates the shear modulus for target models 1 and 3, showing conditions for an exact constant profile, and a profile with some inherent variability, respectively. The Poisson's ratio is set to 0.1 for small-strain applications [56]. For the steel piles, Young's modulus and density are specified as 200 GPa and 7850 kg m⁻³, respectively. The stiffness profiles for each target model are adjusted by applying wₖ, and soil masses are determined through wₘ applied to Eq. (5). A damping ratio of 1% is adopted to represent energy dissipation in the target piles. These piles are considered to be open-ended and fully coring, with the internal soil level matching the external ground level. Additionally, the internal sand's density is assumed to be 2000 kg m⁻³.

In each model, a 10 kN impact load is applied to the pile for a short, arbitrary duration, following which the pile model is allowed to oscillate freely for 50 seconds. The time step for this analysis is set at 0.005 seconds. Each target pile is equipped with five accelerometers positioned along its length (node points in numerical model), starting from the pile's head at 1 m intervals. Data from these five 'sensors' are subsequently utilized in the subspace identification process to determine four natural frequencies of the pile, along with their corresponding mode shapes. These identified parameters serve as the basis for formulating the objective function.
<table><caption>Table 1: Target models</caption><thead><tr><th>Target Model</th><th>w<sub>m</sub></th><th>w<sub>k</sub></th><th>r (m)</th><th>t<sub>0</sub> (m)</th><th>L<sub>1</sub> (m)</th><th>L/D</th><th>Variability<br>in G<sub>0</sub><br>(Y/N)</th></tr></thead><tbody><tr><td>1</td><td>5</td><td>0.5</td><td>0.5</td><td>0.025</td><td>20</td><td>20</td><td>N</td></tr><tr><td>2</td><td>5</td><td>0.5</td><td>2</td><td>0.1</td><td>20</td><td>5</td><td>N</td></tr><tr><td>3</td><td>5</td><td>0.5</td><td>2</td><td>0.1</td><td>20</td><td>5</td><td>Y</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 4: Shear modulus for (a) Target 1 and (b) Target 3. Target 1 represents a constant profile, where Target 3 shows variability in the profile.</figcaption></figure>

As an initial investigation, the objective function was plotted for a number of parameter combinations to examine whether a global minimum existed. Specifically, $w_k$ was set to 0.5 for target model 1, while $L_1$ and $w_m$ were varied. As expected, the objective function, formulated based on the frequencies and mode shapes derived from the five sensors, was minimized for the target parameters, where $w_m = 5$ and $L_1 = 20$ m.
<figure><img src="image_7.png" /><figcaption>Figure 5: Objective function for the parameter pair $L_1$ and $w_m$ in the numerical application for target model 1.</figcaption></figure>

To assess the optimization's performance to estimate $w_m$, $w_k$ and $L_1$ effectively, 100 Monte-Carlo simulations were conducted for each target model. The histograms illustrating the distributions of the estimated parameters are shown in Figs. 6, 7, and 8 for target models 1, 2, and 3, respectively. These histograms represent the probability density functions (PDFs) of the estimated parameters, providing insight into how frequently certain parameter values occur across the simulations. A summary of the results, including the means and standard deviations of the estimated parameters, is presented in Table 2, providing a quantitative assessment of the optimization's performance. From the results, it can be observed that the optimization consistently converges to the expected values for all three parameters. It is noteworthy that in target model 3, where the stiffness has variability, the convergence results in a slightly increased stiffness. This can be explained by the fact that variability in the stiffness generates some bias, causing the actual $w_k$ to shift slightly to the right.
<figure><img src="image_8.png" /><figcaption>Figure 6: Histograms for target 1 showing the distribution of converged values for the three key parameters: (a) $L_1$, (b) $w_m$, and (c) $w_k$. The distributions are obtained from Monte-Carlo simulations and are plotted against the actual values of the corresponding parameters.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Histograms for target 2 showing the distribution of converged values for the three key parameters: (a) $L_1$, (b) $w_m$, and (c) $w_k$. The distributions are obtained from Monte-Carlo simulations and are plotted against the actual values of the corresponding parameters.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 8: Histograms for target 3 showing the distribution of converged values for the three key parameters: (a) $L_1$, (b) $w_m$, and (c) $w_k$. The distributions are obtained from Monte-Carlo simulations and are plotted against the actual values of the corresponding parameters.</figcaption></figure>

<table><caption>Table 2: Mean values and standard deviations of the parameters for target models 1, 2 and 3.</caption><thead><tr><th>Target</th><th>Parameter</th><th>μ</th><th>σ</th><th>Actual Value</th></tr></thead><tbody><tr><td rowspan="3">Target 1</td><td>L<sub>1</sub> (m)</td><td>20.0</td><td>0.017</td><td>20.0</td></tr><tr><td>w<sub>m</sub></td><td>5.0</td><td>0.003</td><td>5.0</td></tr><tr><td>w<sub>k</sub></td><td>0.5</td><td>0.0003</td><td>0.5</td></tr><tr><td rowspan="3">Target 2</td><td>L<sub>1</sub> (m)</td><td>20.0</td><td>0.029</td><td>20.0</td></tr><tr><td>w<sub>m</sub></td><td>5.0</td><td>0.048</td><td>5.0</td></tr><tr><td>w<sub>k</sub></td><td>0.5</td><td>0.001</td><td>0.5</td></tr><tr><td rowspan="3">Target 3</td><td>L<sub>1</sub> (m)</td><td>20.0</td><td>0.023</td><td>20.0</td></tr><tr><td>w<sub>m</sub></td><td>5.0</td><td>0.037</td><td>5.0</td></tr><tr><td>w<sub>k</sub></td><td>0.502</td><td>0.0009</td><td>0.502</td></tr></tbody></table>

## 4. Field Application

Data from a field test reported by Refs. [47, 57] are employed to assess the performance of the optimization approach when applied to real data from impact tests conducted on piles. This assessment utilizes small-strain soil stiffness data, including shear wave velocity measurements and the derived $G_0$ profile. The field test was executed at a quarry in Blessington, southwest of Dublin, Ireland and a general schematic of the pile and soil arrangement is shown in Fig. 1. The experiment involved lateral impact testing of a 7.2 m open-ended steel pile with a diameter of 0.34 m and wall thickness of 0.014 m, embedded in dense, over-consolidated sand. Prior to testing, the soil was excavated leading to an embedment of 4.5m. In the
experiment, the pile is impacted laterally and accelerations are measured using three accelerometers placed
near the top of the pile. For more information on the field testing, refer to Refs. [47, 57].

The quarry at Blessington has undergone extensive characterization in prior research [58] and has been
used to assess the performance of various model, prototype, and full-scale foundation installations, over a
number of years. Comprehensive details on the geotechnical properties of the site can be found in several
studies [59–61]. This field application relies on the site’s small-strain stiffness properties, which were mea-
sured using Multi-Channel Analysis of Surface Waves (MASW) as detailed in Donohue et al. [62]. The
shear wave velocity profile, depicted in Fig. 9a, is used to calculate the small-strain Young’s modulus pro-
file, shown in Fig. 9b. This calculation begins with the small-strain shear modulus ($G_0$) using the formula
$G_0 = \rho v_s^2$, where $v_s$ represents the shear wave velocity. Subsequently, the small-strain Young’s modulus
($E_0$) is calculated using $E_0 = 2G_0(1 + \nu_s)$, where $\rho$ is the soil density (kg m$^{-3}$) and $\nu_s$ is the small-strain
Poisson ratio, assumed here to be 0.1. The site profile data begins at approximately 4.5 m below ground
level due to the excavation process as described above, prior to dynamic testing [47].

<figure><img src="image_11.png" /><figcaption>Figure 9: Small-strain MASW data.</figcaption></figure>

The pile was equipped with three accelerometers along the exposed section of the pile, each set to
record at a frequency of 1000 Hz. The testing protocol involved subjecting the pile to lateral impacts with
a PCB Piezotronics 086D50 model sledgehammer-type modal hammer [63], which features a tip mass of
5.5 kg and recording the resultant impact force time-history. Multiple impacts were performed to assess the
repeatability of the results.
The procedure outlined in Section 2.1 is followed to construct the numerical model for the purpose of model-updating, and Eq. (2) is applied to convert the $G_0$ profile into discrete spring stiffness values $K_n$. The response output data from two impact tests, 1 and 2, are presented in Fig. 10, which shows the acceleration recorded by the three accelerometers located along the pile.

<figure><img src="image_12.png" /><figcaption>Figure 10: Acceleration profiles for impact tests 1 and 2.</figcaption></figure>

The response output data under impact allows for application of different identification methods for modal parameter estimation. To this end, the modal parameters and their variances are identified using the output information-only with the covariance driven subspace identification method [21]. The estimates are obtained for model orders ranging from 20 to 100, 25 time lags to compute the block Hankel matrix and 200 blocks to evaluate its sample covariance. The stabilization diagram of the natural frequency estimates for impacts 1, 2 and 3 is depicted in Fig. 11 alongside the singular value decomposition (SVD) plot, which displays the singular values of the cross power spectral density (CPSD) matrices derived from the acceleration data for impact 2. It can be viewed that among the spurious estimates, two modal alignments at around 12 Hz and 62 Hz can be distinguished, whose frequencies correspond the first and third bending modes for each pile.
<figure><img src="image_13.png" /><figcaption>Figure 11: Stabilization plot to identify excited modes from the field data across impacts 1 to 3 against a power spectral density (PSD) plot from the acceleration data of impact 2.</figcaption></figure>

The function $F(\theta)$ is plotted for the parameter pair ($w_m, L_0$) within the intervals $\{0.5-30\}$ and $\{5-12\}$, respectively, and can be seen in Fig. 12. The objective function is derived from the 1st and 3rd modes of vibration, with mode shapes based on the three accelerometers positioned on the pile. From the objective function, we observe a global minimum around the expected value for the pile length, which is 7.2 m, when $w_k$ is fixed at the mean value of the converged results shown in Table 3. Specifically, in this analysis, $w_k$ was set to the mean value obtained from the Monte Carlo simulations of Impact 1. This global minimum forms a ridge, rather than a well-defined minimum, likely due to the inclusion of only two modes and data from three sensors. This theory is reinforced by the objective function plot for target 1 in Fig. 5, where four modes and data from five sensors are taken into account in the objective function, showing a clearly defined global minimum. This comparison highlights the impact of limited modal information on the optimization landscape.
<figure><img src="image_14.png" /><figcaption>Figure 12: Objective function for the parameter pair $L_0$ and $w_m$ in the field application.</figcaption></figure>

To assess the performance of the optimization to estimate $w_m$, $w_k$ and $L_0$ with real data, Monte Carlo simulations were conducted, similar to the approach described in Section 3. The statistics for $L_0$, $w_m$, and $w_k$ are summarized in Table 3. The analysis reveals that $L_0$ converges to approximately 7.8 m for Impact 1, 7.7 m for Impact 2, and 7.8 m for Impact 3, with standard deviations of 0.48 m, 0.56 m, and 0.55 m, respectively. Similarly, $w_m$ has a mean value of 6.35 with a standard deviation of 2.32 for Impact 1, a mean value of 6.18 with a standard deviation of 3.17 for Impact 2, and a mean value of 6.4 with a standard deviation of 2.65 for Impact 3. For $w_k$, the mean value is 1.73 with a standard deviation of 0.48 for Impact 1, a mean value of 1.59 with a standard deviation of 0.28 for Impact 2, and a mean value of 1.68 with a standard deviation of 0.22 for Impact 3.
<table><caption>Table 3: Mean values and standard deviations of the parameters for Impact 1, 2 and 3.</caption><thead><tr><th>Impact</th><th>Parameter</th><th>μ</th><th>σ</th><th>Actual Value</th></tr></thead><tbody><tr><td rowspan="3">Impact 1</td><td>L<sub>0</sub> (m)</td><td>7.8</td><td>0.48</td><td>7.2</td></tr><tr><td>w<sub>m</sub></td><td>6.35</td><td>2.32</td><td>-</td></tr><tr><td>w<sub>k</sub></td><td>1.73</td><td>0.48</td><td>-</td></tr><tr><td rowspan="3">Impact 2</td><td>L<sub>0</sub> (m)</td><td>7.7</td><td>0.56</td><td>7.2</td></tr><tr><td>w<sub>m</sub></td><td>6.18</td><td>3.17</td><td>-</td></tr><tr><td>w<sub>k</sub></td><td>1.59</td><td>0.28</td><td>-</td></tr><tr><td rowspan="3">Impact 3</td><td>L<sub>0</sub> (m)</td><td>7.8</td><td>0.55</td><td>7.2</td></tr><tr><td>w<sub>m</sub></td><td>6.4</td><td>2.65</td><td>-</td></tr><tr><td>w<sub>k</sub></td><td>1.68</td><td>0.22</td><td>-</td></tr></tbody></table>

The parameter that can be validated is the pile's length, with the exact value being 7.2 m. Although the most likely length as predicted from the procedure is 7.8 m, this is still a good approximation given that the objective function is based on only two modes and three accelerometers positioned close to each other. If more accelerometers were placed along the shaft of the pile and more modes were excited, the objective function would be better defined. This would likely result in the solution converging more closely to the global minimum and the exact value.

Similarly, for $w_m$ and $w_k$, it is observed that the mean values of the converged results are close for impacts 1, 2, and 3. The convergence of all three tests in the same region with similar mean values is significant, as it indicates that different load cases do not interfere with the subspace identification process. More specifically, the identified frequencies and mode shapes for modes 1 and 3 across all impacts yield similar optimization values, reinforcing the robustness and validity of the proposed method.

## 5. Conclusions

This study successfully demonstrates the coupling of model updating, here effectuated via a stochastic evolutionary scheme, the CMA-ES, with subspace identification methods in a robust framework to determine operational mass, stiffness, and embedded length of a number of numerical and experimental piles. The framework adeptly addresses the challenge of estimating critical soil-pile interaction parameters, such as pile length and the associated mobilized mass and stiffness, which are pivotal for accurate foundation behavior modeling. The presented framework combines modal parameter determination through subspace identification with the robustness of the CMA-ES for model updating. While the CMA-ES was chosen due to its practical ability to converge in rugged search landscapes, which may include discontinuities, sharp bends, ridges, noise, local optima, or outliers, any optimization method featuring similar characteristics could be
used. Additionally, CMA-ES features a covariance matrix, reflecting the algorithm’s estimate of the variance
in the search space. A larger covariance matrix indicates a broader exploration, while a smaller covariance
matrix suggests the convergence toward a narrower region, typically near the minimum. Although this fea-
ture was not exploited in the current study, it presents a potential for the uncertainty quantification of the
delivered estimates in the future work. The analysis was conducted in MATLAB on the High-Performance
Computing (HPC) system of the University of Nottingham, specifically utilizing a high-memory compute
node equipped with 2 × 48-core AMD CPUs, 1536 GB of RAM. Under this configuration, the approximate
convergence time for the optimization of every target 1 model was around 12 minutes, with a variance of
1–2 minutes depending on the initial conditions. This approach effectively quantifies uncertainties using
output-only data, enhancing the model’s accuracy, especially for existing piles under operational loads. The
framework’s effectiveness has been validated through numerical model applications and via data reported
from a real full-scale field test, demonstrating its practical applicability and potential to significantly im-
prove soil-structure interaction parameter estimation. Ultimately, this research introduces a framework for
estimating pile parameters using output-only data and highlights the significant potential of statistical up-
dating techniques in advancing the accuracy and reliability of SHM systems in geotechnical and structural
engineering. The approach might be useful for digital twin applications such as for offshore wind turbine
structures.

References

[1] H. J. Burd, C. N. Abadie, B. W. Byrne, G. T. Houlsby, C. M. Martin, R. A. McAdam, R. J. Jardine, A. M. G. Pedro, D. M. Potts, D. M. G. Taborda, L. Zdravković, M. P. Andrade, Application of the pisa design model to monopiles embedded in layered soils, Géotechnique 70 (2020) 1067–1082. doi:10.1680/jgeot.20.PISA.009.
URL https://www.icevirtuallibrary.com/doi/10.1680/jgeot.20.PISA.009

[2] H. J. Burd, D. M. Taborda, L. Zdravković, C. N. Abadie, B. W. Byrne, G. T. Houlsby, K. G. Gavin, D. J. Igoe, R. J. Jardine, C. M. Martin, et al., Pisa design model for monopiles for offshore wind turbines: application to a marine sand, Géotechnique 70 (11) (2020) 1048–1066.

[3] P. Doherty, K. Gavin, Laterally loaded monopile design for offshore wind farms, Proceedings of the institution of civil engineers-energy 165 (1) (2012) 7–17.

[4] B. Gupta, D. Basu, Applicability of timoshenko, euler-bernoulli and rigid beam theories in analysis of laterally loaded monopiles and piles, Géotechnique 68 (9) (2018) 772–785.

[5] E. Winkler, Die Lehre von Elastizität und Festigkeit (on Elasticity and Fixity), Dominicus, 1867.

[6] W. H. Wu, L. J. Prendergast, K. Gavin, An iterative method to infer distributed mass and stiffness profiles for use in reference dynamic beam-winkler models of foundation piles from frequency response functions, Journal of Sound and Vibration 431 (2018) 1–19. doi:10.1016/j.jsv.2018.05.049.

[7] L. J. Prendergast, M. P. Limongelli, N. Ademovic, A. Anžlin, K. Gavin, M. Zanini, Structural health monitoring for performance assessment of bridges under flooding and seismic actions, Structural Engineering International 28 (2018) 296–307. doi:10.1080/10168664.2018.1472534.

URL https://www.tandfonline.com/doi/full/10.1080/10168664.2018.1472534
[8] K. Tatsis, V. Dertimanis, I. Abdallah, E. Chatzi, A substructure approach for fatigue assessment on wind turbine support structures using output-only measurements, Procedia engineering 199 (2017) 1044-1049.

[9] C. R. Farrar, K. Worden, Structural health monitoring: a machine learning perspective, John Wiley & Sons, 2012.

[10] J. Brownjohn, Structural health monitoring of civil infrastructure, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 365 (2007) 589-622. doi:10.1098/rsta.2006.1925.

[11] E. N. Chatzi, C. Papadimitriou, Identification methods for structural health monitoring, Vol. 567, Springer, 2016.

[12] D. J. Ewins, Modal testing: theory, practice and application, John Wiley & Sons, 2009.

[13] S. D. Dhandole, S. V. Modak, A comparative study of methodologies for vibro-acoustic fe model updating of cavities using simulated data, International Journal of Mechanics and Materials in Design 6 (2010) 27-43. doi:10.1007/s10999-010-9116-4
URL http://link.springer.com/10.1007/s10999-010-9116-4

[14] E. Reynders, System identification methods for (operational) modal analysis: Review and comparison, Archives of Computational Methods in Engineering 19 (1) (2012) 51-124.

[15] A. Strauss, S. Kessler, M. P. Limongelli, K. Bergmeister, Existing concrete structures: Life management, testing and structural health monitoring (2020).

[16] B. Peeters, G. de Roeck, Reference-based stochastic subspace identification for output-only modal analysis, Mechanical Systems and Signal Processing 13 (6) (1999) 855-878.

[17] E. Reynders, G. D. Roeck, Reference-based combined deterministic-stochastic subspace identification for experimental and operational modal analysis, Mechanical Systems and Signal Processing 22 (3) (2008) 617-637.

[18] D. Foti, V. Gattulli, F. Potenza, Output-only identification and model updating by dynamic testing in unfavorable conditions of a seismically damaged building, Computer-Aided Civil and Infrastructure Engineering 29 (2014) 659-675. doi:10.1111/mice.12071.

[19] A. Benveniste, L. Mevel, Nonstationary consistency of subspace methods, IEEE Transactions on Automatic Control 52 (6) (2007) 974-984.

[20] E. Reynders, R. Pintelon, G. D. Roeck, Uncertainty bounds on modal parameters obtained from stochastic subspace identification, Mechanical Systems and Signal Processing 22 (2008) 948-969. doi:10.1016/j.ymssp.2007.10.009.

[21] M. Döhler, L. Mevel, Efficient multi-order uncertainty computation for stochastic subspace identification, Mechanical Systems and Signal Processing 38 (2) (2013) 346-366.

[22] P. Mellinger, M. Döhler, L. Mevel, Variance estimation of modal parameters from output-only and input/output subspace-based system identification, Journal of Sound and Vibration 379 (C) (2016) 1-27.

[23] E. P. Reynders, Uncertainty quantification in data-driven stochastic subspace identification, Mechanical Systems and Signal Processing 151 (2021) 107338.

[24] S. Greś, M. Döhler, N.-J. Jacobsen, L. Mevel, Uncertainty quantification of input matrices and transfer function in input/output subspace system identification, Mechanical Systems and Signal Processing 167 (2022) 108581.

[25] M. Friswell, J. E. Mottershead, Finite element model updating in structural dynamics, Vol. 38, Springer Science & Business Media, 1995.

[26] M. I. Friswell, D. J. Inman, D. F. Pilkey, Direct updating of damping and stiffness matrices, AIAA Journal 36 (1998) 491-493. doi:10.2514/2.396.

[27] S. V. Modak, Direct matrix updating of vibroacoustic finite element models using modal test data, AIAA Journal 52 (2014) 1386-1392. doi:10.2514/1.J052558
URL https://arc.aiaa.org/doi/10.2514/1.J052558

[28] S. Sehgal, H. Kumar, Structural dynamic model updating techniques: A state of the art review, Archives of Computational Methods in Engineering 23 (2016) 515-533. doi:10.1007/s11831-015-9150-3.
URL http://link.springer.com/10.1007/s11831-015-9150-3

[29] N. F. Alkayem, M. Cao, Y. Zhang, M. Bayat, Z. Su, Structural damage detection using finite element model updating with evolutionary algorithms: a survey (7 2018). doi:10.1007/s00521-017-3284-1.

[30] D. Nehete, S. Modak, K. Gupta, Structural fe model updating of cavity systems incorporating vibro-acoustic coupling, Mechanical Systems and Signal Processing 50-51 (2015) 362-379. doi:10.1016/j.ymssp.2014.05.028.
URL https://linkinghub.elsevier.com/retrieve/pii/S0888327014002040

[31] F. Sarvi, S. Shojaee, P. Torkzadeh, Damage identification of trusses by finite element model updating using an enhanced levenberg-marquardt algorithm (2014).

[32] S. Modak, T. Kundra, B. Nakra, Model updating using constrained optimization, Mechanics Research Communications 27 (2000) 543-551. doi:10.1016/S0093-6413(00)00128-2.

[33] S. Pradhan, S. Modak, A method for damping matrix identification using frequency response data, Mechanical Systems and Signal Processing 33 (2012) 69-82. doi:10.1016/j.ymssp.2012.07.002.
URL https://linkinghub.elsevier.com/retrieve/pii/S0888327012002658

[34] A. Ioakim, L. Prendergast, A two-stage method to estimate the embedded length of foundation piles using frf-based model updating, Mechanical Systems and Signal Processing 219 (2024) 111603. doi:10.1016/j.ymssp.2024.111603.

[35] H. Sohn, K. H. Law, A bayesian probabilistic approach for structure damage detection, Earthquake engineering & structural dynamics 26 (12) (1997) 1259-1281.

[36] X. Jiang, S. Mahadevan, Bayesian probabilistic inference for nonparametric damage detection of structures, Journal of engineering mechanics 134 (10) (2008) 820-831.

[37] A. H. Alavi, H. Hasni, N. Lajnef, K. Chatti, F. Faridazar, An intelligent structural damage detection approach based on self-powered wireless sensor data, Automation in Construction 62 (2016) 24-44. doi:10.1016/j.autcon.2015.10.001.
URL https://linkinghub.elsevier.com/retrieve/pii/S092658051500206X

[38] T. Marwala, Finite element model updating using computational intelligence techniques: applications to structural dynamics, Springer Science & Business Media, 2010.

[39] S. Greś, M. Döhler, L. Mevel, Statistical model-based optimization for damage extent quantification, Mechanical Systems and Signal Processing 160 (2021) 107894.

[40] S. Greś, R. Riva, C. Y. Süleyman, P. Andersen, M. Luczak, Uncertainty quantification of modal parameter estimates obtained from subspace identification: An experimental validation on a laboratory test of a large-scale wind turbine blade, Engineering Structures 256 (2022) 114001.

[41] J. D. Holt, S. Chen, R. A. Douglas, Determining lengths of installed timber piles by dispersive wave propagation, Transportation research record (1994).

[42] A. A. Hanifah, A theoretical evaluation of guided waves in deep foundations, Northwestern University, 1999.

[43] H.-c. Chao, An experimental model for non-destructive evaluation on pile foundations using guided wave approach, Ph.D. thesis, Northwestern University (2002).

[44] H. Wang, Theoretical evaluation of embedded plate-like and solid cylindrical concrete structures with guided waves, Northwestern University, 2004.

[45] V. Samu, M. Guddati, Nondestructive method for length estimation of pile foundations through effective dispersion analysis of reflections, Journal of Nondestructive Evaluation 38 (6 2019). doi:10.1007/s10921-019-0583-8.

[46] S. Cui, H. Liu, P. Maghoul, Pile length estimation based on guided wave theory for unknown foundations, in: The 75th Annual Canadian Geotechnical Society Conference (GeoCalgary 2022), 2022, pp. 1-6.

[47] L. J. Prendergast, K. Gavin, A comparison of initial stiffness formulations for small-strain soil-pile dynamic winkler modelling, Soil Dynamics and Earthquake Engineering 81 (2016) 27-41.

[48] D. Kallehave, C. L. Thilsted, M. Liingaard, Modification of the api py formulation of initial stiffness of sand, in: SUT
Offshore Site Investigation and Geotechnics, SUT, 2012, pp. SUT-OSIG.

[49] S. Timoshenko, History of strength of materials: with a brief account of the history of theory of elasticity and theory of structures, Courier Corporation, 1983.

[50] A. B. Vesić, Bending of beams resting on isotropic elastic solid, Journal of the Engineering Mechanics Division 87 (1961) 35-53. doi:10.1061/JMCEA3.0000212.

[51] M. Döhler, L. Mevel, Fast multi-order computation of system matrices in subspace-based system identification, Control Engineering Practice 20 (9) (2012) 882–894.

[52] D. Bauer, M. Deistler, W. Scherrer, Consistency and asymptotic normality of some subspace algorithms for systems without observed inputs, Automatica 35 (7) (1999) 1243–1254.

[53] E. Reynders, J. Houbrechts, G. De Roeck, Fully automated (operational) modal analysis, Mechanical Systems and Signal Processing 29 (2012) 228–250.

[54] N. Hansen, A. Ostermeier, Completely derandomized self-adaptation in evolution strategies, Evolutionary Computation 9 (2) (2001) 159–195.

[55] J. Tott-Buswell, L. J. Prendergast, K. Gavin, A cpt-based multi-spring model for lateral monopile analysis under sls conditions in sand, Ocean Engineering 293 (2024) 116642.

[56] L. J. Prendergast, D. Hester, K. Gavin, J. O'Sullivan, An investigation of the changes in the natural frequency of a pile affected by scour, Journal of sound and vibration 332 (25) (2013) 6685-6702.

[57] L. Prendergast, W. Wu, K. Gavin, Experimental application of frf-based model updating approach to estimate soil mass and stiffness mobilised under pile impact tests, Soil Dynamics and Earthquake Engineering 123 (2019) 1-15.

[58] P. Doherty, K. Gavin, Laterally loaded monopole design for offshore wind farms, Proceedings of Institution of Civil Engineers: Energy 165 (2012) 7–17. doi:10.1680/ener.11.00003.

[59] L. J. Prendergast, D. Hester, K. Gavin, J. J. O'Sullivan, An investigation of the changes in the natural frequency of a pile affected by scour, Journal of Sound and Vibration 332 (2013) 6685-6702. doi:10.1016/j.jsp.2013.08.020.

[60] K. Gavin, B. Lehane, Base load – displacement response of piles in sand, Canadian Geotechnical Journal 44 (2007) 1053-1063. doi:10.1139/T07-048.

[61] K. Gavin, A. Adekunte, B. O'Kelly, A field investigation of vertical footing response on sand, Proceedings of the Institution of Civil Engineers - Geotechnical Engineering 162 (2009) 257-267. doi:10.1680/geng.2009.162.5.257.

[62] S. Donohue, M. Long, K. Gavin, P. O'Connor, Shear wave stiffness of irish glacial till, in: Proceedings of the 2nd International Conference on Site Characterization, 2004, pp. 459–466.

[63] P. Piezotronics, Impulse hammer 086d50. model 086d50 prod specif (2018).