Assessing the small-strain soil stiffness for offshore wind turbines based on in situ seismic measurements

The fundamental natural frequency as measured on installed offshore wind turbines is significantly higher than its designed value, and it is expected that the explanation for this can be found in the currently adopted modeling of soil-structure interaction. The small-strain soil stiffness is an important design parameter, as it has a defining influence on the first natural frequency of these structures. In this contribution, in situ seismic measurements are used to derive the small-strain shear modulus of soil as input for 3D soil-structure interaction models to assess the initial soil stiffness at small strains for offshore wind turbine foundations. A linear elastic finite element model of a half-space of solids attached to a pile is used to derive an equivalent first mode shape of the foundation. The second model extends the first one by introducing contact elements between pile and soil, to take possible slip and gap-forming into account. The deflections derived with the 3D models are smaller than those derived with the p- y curve design code. This higher stiffness is in line with the higher measured natural frequencies. Finally a method is suggested to translate the response of 3D models into 1D engineering models of a beam laterally supported by uncoupled distributed springs.


Introduction
In order to become an economically competitive energy source, the levelized cost of electricity generated by offshore wind power needs to be lowered. Less conservatism in design of the support structure (SStr) of offshore wind turbines (OWT) can contribute to this by decreasing the capital cost and/or increasing the energy yield. The latter is achieved by a longer certified lifetime of already installed OWTs. Both can be attained by increased knowledge of soil-structure interaction (SSI): in depth understanding of both the stiffness and damping occurring in the interaction of the monopile and the soil.
Over the recent years, measured natural frequencies of installed OWTs have been found to be higher than designed for. Apart from being a waste of expensive steel -often larger diameter foundation piles are employed for reaching the desired fundamental natural frequency -also more fatigue damage is accumulated than anticipated, because of the higher amount of cycles during its lifetime. The industry is turning towards SSI for finding the cause of this discrepancy in frequency between design and real life. A higher measured frequency may suggest an under-estimation of the stiffness of the soil against the pile motion.
The current engineering approach which is prescribed in the standards [1] is called the "p-y curve" method. In this approach, the monopile is laterally supported by uncoupled distributed non-linear springs. The displacement 'y'-dependent magnitude of these springs is based on empirical relations found for slender, small-diameter piles in both clay [2] and sand [3]. However, the ratio of embedded length L over diameter D that is currently employed in the industry, is expected to invoke a fundamentally different soil reaction than for which the p-y curves were originally calibrated. The typical monopile bending can be characterized by being more rigid than the flexible piles for which the p-y curves are to be used. Different soil-reaction mechanisms are associated with the bending behavior of monopiles than with the bending of slender flexible piles. Apart from this, it is expected to be wrong to assume that the magnitude of the initial stiffness is independent of the diameter of the pile, and that this stiffness would increase linearly with depth (as is done in the p-y curve method) [4].
The initial stiffness has a large effect on the natural frequency of the OWT: the soil responds linearly for most of the endured vibrational amplitudes of the foundation, and this small-strain stiffness defines the modal properties of the SStr. Therefore, this initial part of the reaction, often called initial subgrade modulus k S,0 or E* py [N/m 2 ], has been studied by several researchers [5], [6], [7]. The dimension of k s,0 can be understood as a discretized spring stiffness over a unit vertical length of the pile, so it is rather [N/m/m]. In the p-y curve formulation k s,0 is calculated by multiplying k, the modulus of subgrade reaction [N/m 3 ] by depth (z). k is determined with the angle of internal friction of sand φ.
However, seemingly contradictory to what is mentioned above, these referred studies find that the p-y curve method actually over-predicts the initial stiffness of the soil for larger depths. If this were generally the case, it would not explain the higher natural frequencies that are measured. The researchers suggest that a more realistic -taking 3D, global effects into account -variation of the stiffness with depth might rather follow a power law form with a power exponent smaller than 1 [4].
Geotechnical tests like the cone penetration tests (CPT) and borehole samples for laboratory testing are the most applied soil measurements in the wind industry. A balanced approach with geophysical (seismic) measurements, with a focus on the low frequency response, would be beneficial. Whereas the CPT measures local resistance and friction along the shaft, a recorded soil wave carries information of the soil characteristics along its entire path, yielding more global parameters of a certain location. This global parameterization of the soil can best be used to predict the response of the foundation. The in situ stiffness is determined by assessing the wave velocities in the profile. In this way, all local properties are incorporated, like for instance oedometric stiffening and saturation effects. This can be considered an advantage over most laboratory tests in which these in situ characteristics are partially lost.
We aim to assess this initial stiffness aspect. An alternative approach is suggested, in which a 3D model is used to calibrate the initial part of the p-y curves which are fit for fast simulation design purposes. This initial reaction, like is done in the p-y curve method, is assumed to be linear elastic. The 3D model incorporates global effects and interface friction. The novelty of the approach lies in using in situ seismic measurements to extract the small-strain stiffness of the design location. The material properties identified with the seismic measurements are input to the 3D model. Finally, a translation method is suggested for finding the 1D equivalent stiffness (i.e., the initial branch of the py curve) of the 3D response.

Seismic in situ data
An in situ soil measurement campaign was carried out at a near shore wind farm. This campaign comprised a geotechnical part including a CPT at each OWT location and 10 boreholes for lab test sampling, followed by a geophysical campaign. At 6 locations around the farm, seismic measurements including SCPTs and a newly developed test designed to capture the low-frequency dynamic response of the soil: the Low Frequency Cone Penetration Tests (LFCPT) were performed. This paper will discuss a frequency independent analysis of the SCPT data.
Here we focus on the seismic data of a location closest to one of the design locations (deepest & softest soil). A depth of 25m was reached at this location, measuring each meter with a dual-phone cone with an interval distance of 0.5m. A hydraulic shear wave hammer was used as excitation device at mudline. Stacking responses of multiple hits for each depth rendered clear shear-wave patterns. To enable automatic picking, the wave arrival time was defined as the maximum peak.
The difference in arrival time between two adjacent receivers -the interval time -is obtained by simply subtracting the arrival time of the upper geophone from the lower geophone. For example, for the third layer in Figure 1, this is Δt = t 32 -t 31 . A minimization problem was set up to find the shear wave velocities. Assuming a stratification with 1m layer thickness, the successive layer shear wave velocities can be computed by minimizing an objective function for the observed arrival time. The objective function incorporates the effect of wave refraction through Snell's law. Confidence in finding the global minimum can be reached by visual inspection of the objective function for the first layers, and is aided by a realistic initial guess of the solution.
A schematic view of the solution method is shown in Figure 1, in which an example is given for the third layer. Geometric relations with their unknowns form the constraining equations of the optimization problem. The amount of equations and unknowns equals 4n+1, where 'n' denotes the layer number. Figure 1: Assumed wave ray scheme to find wave velocities. An example is given for layer number 3. φ inc is the incoming wave angle, φ refr is the angle of the refracted wave at the assumed layer interface The determined shear wave velocities are displayed in Figure 2. Next to the velocity profile, the laboratory soil classification is given, which is based on samples retrieved from a borehole at the same location. It can be seen that the velocity increases with depth, which is expected. There are some outliers, and the weaker soils between approximately 15m and 23m depth are reflected in the measured velocities.
The in situ soil densities, ρ, were determined using the cone tip and frictional resistance measured with the SCPT. According to the type of soil, a Poisson's ratio, ν, was estimated: 0.3 for sand and 0.45 for cohesive material. The Young's modulus, E, was calculated according to in which V s is the shear wave velocity, and G the shear modulus.

Modeling approach
Two FE models were developed using ANSYS to try and capture the global (non-local) stiffness effects of the interaction between pile and soil. The models are similar in terms of the soil and the pile, but differ in the way the interface between these two media is modeled. For both models, the pile was modeled with Shell elements, and the soil was modeled using Solid elements. The dimensions of the pile and halfspace are equal for both models. The pile has a diameter of 5m, embedded length of 32m and a wall thickness of 60mm. The soil halfspace was given a vertical dimension of 50m, and a radius of 100m. An impression of the model is given in Figure 3. The models are described in the next 2 paragraphs. In paragraph 4.2 a method is proposed to find a 1D equivalent lateral stiffness for the one found with a 3D model.

Linear elastic FE model
The first model is linear elastic. The soil elements in and outside the pile, are attached to the pile elements, creating the possibility of high shear forces on the interface, and tensile forces between the soil and the pile. However, as the focus is laid upon small-amplitude vibrations, such kinematicphysically incorrect -constraints might not contaminate the results too much. A unit horizontal force of 1 MN and bending moment of 1 MNm was applied to the top of the pile at 5m above mudline. The soil stratification was given a 1m discretization and each layer was assigned the material properties derived from the seismic and CPT data as described in section 2. Because of the limited depth of the seismic data, the deepest 7m along the pile (from -25 to -32m) was assumed to be one layer, with the same properties as the layer above (-24 to -25m). Geotechnical data indicates that indeed a quite uniform sand layer is present at this depth. However, the assumption of constant Young's modulus can be conservative, as it could be higher due to a larger effective overburden pressure.

Elastic FE model including contact condition
The second model is almost identical to the linear elastic model, but includes the use of Contact Elements (CE) at the interface of pile and soil. Gaps between the soil and the pile can now occur, and a maximum friction coefficient, μ can be defined. If the maximum friction force is exceeded, sliding will occur.
The interface friction angle δ was taken to be 2/3 of the internal angle of friction of the soil. Then μ was calculated as tan(δ). Averaging this over the depth results in μ=0.5. Three cases were analyzed: firstly a homogeneous soil was assumed with a friction coefficient of 0.5. The soil properties of the sand layer between -13m and -14m (E=179Mpa, ν=0.3 and ρ=1800kg/m) 3 were given to the entire soil stratum for this homogeneous case, as these values are close to the average soil properties over the embedded depth. Then the same friction coefficient was taken, but with the 'real' stratified soil. Thirdly a stratified case with friction coefficient μ=0.4 was calculated for comparison purposes.
Summing up the three cases presented are: • A homogeneous case in which an average density, Young's modulus and Poisson's ratio were used, and μ=0.5 • A stratified ('realistic') case with μ=0.5 • A stratified case with μ=0.4 With the limited forcing that is applied, the soil 'sticks' to the pile nearly everywhere, except for some sliding in and outside the pile near mudline. This generally holds for the three cases.

Results
The deflection shapes and lateral soil reactions from the 3D models are presented in the next paragraph. A method is suggested for finding a local, uncoupled equivalent spring stiffness distribution in paragraph 4.2. Figure 4 displays the deflection shapes and zoom of the pile-tip displacements for the linear elastic model, the design (p-y) code model and the 3 CE cases. The grey line reflects the shape derived when applying the design code [1], and applying double the load which was used for the halfspace models (a horizontal force of 2MN and moment of 2MNm), as this 1D model does not consider half of the symmetric problem. The black line, 'Linear FE model', is the shape of the linear elastic model described in paragraph 3.1. The other 3 shapes are the cases described in paragraph 3.2 including interface friction (variation) and an average homogeneous soil case.

Deflections and horizontal reactions derived from 3D model
It can be seen that the p-y curve approach seems to be conservative in estimating the displacements at mudline when compared to the 3D FE models. The main reasons for this are thought to be twofold; the input for the FE models was derived with seismic measurements, which is an appropriate way for determining the small strain shear modulus. Secondly, these FE models incorporate 3D phenomena for soil resistance to the large diameter piles. These phenomena can for instance include pressure redistribution in the soil due to the Poisson's effect and arching. The pile reaction modeled with the 1D p-y method behaves more flexible, and especially the upper 10m (30% of the depth) deflections are significantly larger than the global models predict, resulting in up to 50% larger deflections at mudline. From the stiffness distribution displayed in Figure 6, we can see that the grey line from the p-y curve approach would indeed follow a linearly increasing initial stiffness with depth, despite the soft soil layer between -14m and -23m. This creates high stiffness at greater depth, explaining the rather small 'toe-kick' (lateral displacement of pile tip) for this design approach.
The deflection shape of the linear FE model has a smaller mudline displacement than the design code, and its toe-kick is the greatest of all models. As previously mentioned, because of the attachment of the soil Solid elements to the Shells of the pile, this model allows higher interface friction forces than the models including a slip possibility.
Further, we see that the shapes of the two cases of contact element FE models with friction coefficients μ of 0.4 (indicated with '+' points) and 0.5 (the '-·' line) are quite similar. The point where these curves cross "zero-displacement" is shifted upwards with respect to the black solid line of the linear model. Because of smaller soil pressures near mudline, the models with contact elements allow for some slip in the upper region, which might explain the larger mudline deflections.
Finally, the dashed ('--') line reflects the deflection shape when assuming a homogeneous soil with the averaged properties of the true stratified soils, and a friction coefficient of μ=0.5. The shape is quite different which indicates that stratification has a large influence on the response. This can be expected: an averaged stiffness gives significantly smaller displacements at both mudline as at the pile-tip. The presence of a stiffer soil layer (because of averaging) at a shallow depth, and possibly also around the point of rotation, has a large influence on the displacements at mudline and consequently also at the pile-tip. The zero-displacement crossing is identical to the stratified contact models.  When summing the horizontal nodal forces along the shells of the piles (per 1 m vertical length), the lateral force of the soil, p [N/m], on the pile are derived as displayed in Figure 5. When subsequently dividing these forces p, over the displacements y of Figure 4, we get a type of local initial stiffness, k s,0 [N/m 2 ] as shown in Figure 6.  Figure 6 it can be concluded that the initial stiffness obtained from the global 3D models is not the equivalent initial stiffness that is to be used for the 1D engineering models. Due to the singularity effects when dividing the force over (near-) zero displacements, none-physical high stiffnesses are found, and also negative stiffness can occur above the zero-crossing. The latter might be caused by a global effect: soil laterally displaced below the zero-crossing, can increase the lateral stress in the soil above this point due to the Poisson's effect. A similar high stiffness due to the zerocrossing was also found by Sørensen et al [7]. In the next paragraph, a method is suggested to translate the results found with the more realistic 3D global models in a 1D engineering model.

Translation of 3D response to 1D engineering model
In order to use the results of (more advanced) global models for simpler engineering models, one is faced with the challenge of capturing the global effects of a 3D model into a 1D model. A simple division of the lateral force, p [N/m] over y [m] does not give useable 1D uncoupled spring stiffnesses, as shown in paragraph 4.1.
To find the equivalent p-y curve initial stiffness, we propose to optimize the stiffness profile in the engineering model to match the deflection shape obtained from a more advanced model. To this end, a 1D FE model was used to model a Timoshenko beam, laterally supported by uncoupled springs. Naturally the same pile properties and pile-head forcing were applied.
The following objective function was defined: ) ) ( min( in which n is the node number of the beam, and N the amount of nodes, u b is the horizontal displacement of the Timoshenko beam model, and u t the target displacement from one of the 3D models. As fitting a deflection shape with a certain stiffness profile does not have a unique solution, this optimization is performed in a stepwise approach in order to steer the solution to be physically acceptable and within the bounds of the classic engineering beam-on-Winkler foundation model, which does not allow negative stiffness. To this end, the first step is to assume a constant stiffness distribution with depth, and finding the initial stiffness which gives the best fit to the target deflection. After that, a linear increasing initial stiffness with depth is assumed. The stiffness distribution with depth is assumed to be a polynomial, of which the order is increased with each step. The four distributions considered in this contribution are given by: (4)(5) in which equations 4-2 -4-5 are respectively the constant, linear, quadratic and cubic stiffness distributions with depth.
Since the constant distribution is a one dimensional minimization problem (only K 11 in equation 4-2 is varied), the solution space can be plotted and it can be visually verified that the global minimum has been found. The resulting deflection and stiffness distribution can be seen in Figure 7, in which also the deflection and initial stiffness distribution is given which is obtained by following the design code. It is clear that a good fit cannot be obtained when assuming constant stiffness with depth. The right panel of Figure 8 also includes the ´local' stiffness of the 3D model, which, as previously discussed in section 4.1, cannot be used for engineering purposes.   Figure 7. Note that K 11 ≠K 21 ≠K 31 nor K 22 ≠K 32 . The resulting deflection fit is displayed in Figure 9 and the associated linear stiffness distribution with depth is displayed in Figure 10. In the same two respective plots, also the deflection shapes and stiffness distributions of the quadratic (equation 4-4) and cubic (equation [4][5] polynomial assumptions are shown. The quality of the fit clearly increases with the order of the Although the fit is improved with higher order of the polynomial, we can question whether a perfect fit can be obtained when trying to mimic the deflection from a global model with only lateral springs. It seems challenging to match the same bending in the pile without the use of, for instance, rotational springs or varying the Young's modulus of the pile along its length. The latter was checked, and the quality of the fit was much enhanced. However, varying the pile's Young's modulus along its length, is not realistic, and cannot be used for engineering purposes.

Conclusions and Discussion
When using the in situ measured small-strain shear modulus in a global 3D model, a first run of analyses indicated smaller mudline displacements than the 'p-y curve' approach dictated by the industry design code [1]. This result may be interpreted as being according to expectation. Although researchers' opinions on this subject deviate [8], the presented method of using in situ determined small-strain properties in a 3D global model was expected to generate higher stiffness. This would explain the higher measured natural frequencies than designed for (see section 1).
A method is suggested for finding a 1D equivalent stiffness distribution of the obtained global responses. In this optimization method, the order of the polynomial of the stiffness distribution with depth was stepwise increased. The quality of the fit improves with each order of the polynomial, and indeed other researchers have gone up to 4 th order fits of lateral soil resistance [9]. However, incorporating only lateral springs (apart from for instance also including rotational springs) seems to limit the resemblance of the 1D model. We should also realise that the solution domain is bounded by the assumption that a polynomial form can be used to describe the stiffness with depth. Some researchers [4] find lower initial stiffness than the p-y curve approach for greater depths, and higher stiffness for shallow depths. The location of the stiff layers in the depth profile is expected to have a high influence on the natural frequency. This is confirmed by the different deflection shape when considering a homogeneous soil, as can be seen in Figure 4. Therefore, we need to adopt a different The presented work is only a part of the investigations that need to be performed. Further research is needed to be able to draw conclusions on the assessment of the small-strain soil stiffness for offshore wind turbines.